]>
Commit | Line | Data |
---|---|---|
0879e280 | 1 | /************************************************************************** |
2 | * Author: Panos Christakoglou. * | |
3 | * Contributors are mentioned in the code where appropriate. * | |
4 | * * | |
5 | * Permission to use, copy, modify and distribute this software and its * | |
6 | * documentation strictly for non-commercial purposes is hereby granted * | |
7 | * without fee, provided that the above copyright notice appears in all * | |
8 | * copies and that both the copyright notice and this permission notice * | |
9 | * appear in the supporting documentation. The authors make no claims * | |
10 | * about the suitability of this software for any purpose. It is * | |
11 | * provided "as is" without express or implied warranty. * | |
12 | **************************************************************************/ | |
13 | ||
14 | /* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */ | |
15 | ||
16 | //----------------------------------------------------------------- | |
17 | // Balance Function class | |
18 | // This is the class to deal with the Balance Function wrt Psi analysis | |
19 | // Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch | |
20 | //----------------------------------------------------------------- | |
21 | ||
22 | ||
23 | //ROOT | |
24 | #include <Riostream.h> | |
621329de | 25 | #include <TCanvas.h> |
0879e280 | 26 | #include <TMath.h> |
27 | #include <TAxis.h> | |
28 | #include <TH2D.h> | |
29 | #include <TH3D.h> | |
30 | #include <TLorentzVector.h> | |
31 | #include <TObjArray.h> | |
32 | #include <TGraphErrors.h> | |
33 | #include <TString.h> | |
34 | ||
35 | #include "AliVParticle.h" | |
36 | #include "AliMCParticle.h" | |
37 | #include "AliESDtrack.h" | |
38 | #include "AliAODTrack.h" | |
9afe3098 | 39 | #include "AliTHn.h" |
35aff0f3 | 40 | #include "AliAnalysisTaskTriggeredBF.h" |
0879e280 | 41 | |
42 | #include "AliBalancePsi.h" | |
c2aad3ae | 43 | using std::cout; |
44 | using std::endl; | |
45 | using std::cerr; | |
0879e280 | 46 | |
47 | ClassImp(AliBalancePsi) | |
48 | ||
49 | //____________________________________________________________________// | |
50 | AliBalancePsi::AliBalancePsi() : | |
51 | TObject(), | |
9fd4b54e | 52 | fShuffle(kFALSE), |
0879e280 | 53 | fAnalysisLevel("ESD"), |
54 | fAnalyzedEvents(0) , | |
55 | fCentralityId(0) , | |
56 | fCentStart(0.), | |
57 | fCentStop(0.), | |
9afe3098 | 58 | fHistP(0), |
59 | fHistN(0), | |
60 | fHistPN(0), | |
61 | fHistNP(0), | |
62 | fHistPP(0), | |
63 | fHistNN(0), | |
73609a48 | 64 | fHistHBTbefore(0), |
65 | fHistHBTafter(0), | |
66 | fHistConversionbefore(0), | |
67 | fHistConversionafter(0), | |
c8ad9635 | 68 | fHistPsiMinusPhi(0), |
f0754617 | 69 | fHistResonancesBefore(0), |
70 | fHistResonancesRho(0), | |
71 | fHistResonancesK0(0), | |
72 | fHistResonancesLambda(0), | |
5a861dc6 | 73 | fHistQbefore(0), |
74 | fHistQafter(0), | |
73609a48 | 75 | fPsiInterval(15.), |
7c04a4d2 | 76 | fDeltaEtaMax(2.0), |
5a861dc6 | 77 | fResonancesCut(kFALSE), |
78 | fHBTCut(kFALSE), | |
d9eb282c | 79 | fHBTCutValue(0.02), |
5a861dc6 | 80 | fConversionCut(kFALSE), |
a9737921 | 81 | fInvMassCutConversion(0.04), |
5a861dc6 | 82 | fQCut(kFALSE), |
6fa567bd | 83 | fDeltaPtMin(0.0), |
683cbfb5 | 84 | fVertexBinning(kFALSE), |
f0fb4ac1 | 85 | fEventClass("EventPlane"){ |
0879e280 | 86 | // Default constructor |
0879e280 | 87 | } |
88 | ||
0879e280 | 89 | //____________________________________________________________________// |
90 | AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance): | |
9fd4b54e | 91 | TObject(balance), fShuffle(balance.fShuffle), |
0879e280 | 92 | fAnalysisLevel(balance.fAnalysisLevel), |
93 | fAnalyzedEvents(balance.fAnalyzedEvents), | |
94 | fCentralityId(balance.fCentralityId), | |
95 | fCentStart(balance.fCentStart), | |
96 | fCentStop(balance.fCentStop), | |
9afe3098 | 97 | fHistP(balance.fHistP), |
98 | fHistN(balance.fHistN), | |
99 | fHistPN(balance.fHistPN), | |
100 | fHistNP(balance.fHistNP), | |
101 | fHistPP(balance.fHistPP), | |
102 | fHistNN(balance.fHistNN), | |
73609a48 | 103 | fHistHBTbefore(balance.fHistHBTbefore), |
104 | fHistHBTafter(balance.fHistHBTafter), | |
105 | fHistConversionbefore(balance.fHistConversionbefore), | |
106 | fHistConversionafter(balance.fHistConversionafter), | |
c8ad9635 | 107 | fHistPsiMinusPhi(balance.fHistPsiMinusPhi), |
f0754617 | 108 | fHistResonancesBefore(balance.fHistResonancesBefore), |
109 | fHistResonancesRho(balance.fHistResonancesRho), | |
110 | fHistResonancesK0(balance.fHistResonancesK0), | |
111 | fHistResonancesLambda(balance.fHistResonancesLambda), | |
5a861dc6 | 112 | fHistQbefore(balance.fHistQbefore), |
113 | fHistQafter(balance.fHistQafter), | |
73609a48 | 114 | fPsiInterval(balance.fPsiInterval), |
7c04a4d2 | 115 | fDeltaEtaMax(balance.fDeltaEtaMax), |
2922bbe3 | 116 | fResonancesCut(balance.fResonancesCut), |
73609a48 | 117 | fHBTCut(balance.fHBTCut), |
d9eb282c | 118 | fHBTCutValue(balance.fHBTCutValue), |
f0fb4ac1 | 119 | fConversionCut(balance.fConversionCut), |
a9737921 | 120 | fInvMassCutConversion(balance.fInvMassCutConversion), |
5a861dc6 | 121 | fQCut(balance.fQCut), |
6fa567bd | 122 | fDeltaPtMin(balance.fDeltaPtMin), |
683cbfb5 | 123 | fVertexBinning(balance.fVertexBinning), |
f0fb4ac1 | 124 | fEventClass("EventPlane"){ |
0879e280 | 125 | //copy constructor |
0879e280 | 126 | } |
127 | ||
128 | //____________________________________________________________________// | |
129 | AliBalancePsi::~AliBalancePsi() { | |
130 | // Destructor | |
9afe3098 | 131 | delete fHistP; |
132 | delete fHistN; | |
133 | delete fHistPN; | |
134 | delete fHistNP; | |
135 | delete fHistPP; | |
136 | delete fHistNN; | |
73609a48 | 137 | |
138 | delete fHistHBTbefore; | |
139 | delete fHistHBTafter; | |
140 | delete fHistConversionbefore; | |
141 | delete fHistConversionafter; | |
c8ad9635 | 142 | delete fHistPsiMinusPhi; |
f0754617 | 143 | delete fHistResonancesBefore; |
144 | delete fHistResonancesRho; | |
145 | delete fHistResonancesK0; | |
146 | delete fHistResonancesLambda; | |
5a861dc6 | 147 | delete fHistQbefore; |
148 | delete fHistQafter; | |
f0fb4ac1 | 149 | |
0879e280 | 150 | } |
151 | ||
152 | //____________________________________________________________________// | |
9afe3098 | 153 | void AliBalancePsi::InitHistograms() { |
154 | // single particle histograms | |
211b716d | 155 | |
156 | // global switch disabling the reference | |
157 | // (to avoid "Replacing existing TH1" if several wagons are created in train) | |
158 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
159 | TH1::AddDirectory(kFALSE); | |
160 | ||
9afe3098 | 161 | Int_t anaSteps = 1; // analysis steps |
9fd4b54e | 162 | Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables |
163 | Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables | |
164 | TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables | |
9afe3098 | 165 | |
166 | // two particle histograms | |
9fd4b54e | 167 | Int_t iBinPair[kTrackVariablesPair]; // binning for track variables |
168 | Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables | |
169 | TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables | |
f0fb4ac1 | 170 | /********************************************************** |
171 | ||
172 | ======> Modification: Change Event Classification Scheme | |
173 | ||
174 | ---> fEventClass == "EventPlane" | |
175 | ||
176 | Default operation with Event Plane | |
177 | ||
178 | ---> fEventClass == "Multiplicity" | |
179 | ||
180 | Work with kTPCITStracklet multiplicity (from GetReferenceMultiplicity) | |
181 | ||
182 | ---> fEventClass == "Centrality" | |
183 | ||
184 | Work with Centrality Bins | |
185 | ||
186 | ***********************************************************/ | |
187 | ||
188 | //--- Multiplicity Bins ------------------------------------ | |
a9f20288 | 189 | const Int_t kMultBins = 10; |
f0fb4ac1 | 190 | //A first rough attempt at four bins |
a9f20288 | 191 | Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80,100,100000}; |
f0fb4ac1 | 192 | //---------------------------------------------------------- |
193 | ||
194 | //--- Centrality Bins -------------------------------------- | |
683cbfb5 | 195 | const Int_t kNCentralityBins = 9; |
a9f20288 | 196 | const Int_t kNCentralityBinsVertex = 15; |
683cbfb5 | 197 | Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; |
a9f20288 | 198 | Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; |
f0fb4ac1 | 199 | //---------------------------------------------------------- |
200 | ||
201 | //--- Event Plane Bins ------------------------------------- | |
202 | //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest) | |
203 | const Int_t kNPsi2Bins = 4; | |
204 | Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5}; | |
205 | //---------------------------------------------------------- | |
206 | ||
207 | //Depending on fEventClass Variable, do one thing or the other... | |
208 | if(fEventClass == "Multiplicity"){ | |
209 | iBinSingle[0] = kMultBins; | |
210 | dBinsSingle[0] = kMultBinLimits; | |
211 | axisTitleSingle[0] = "kTPCITStracklet multiplicity"; | |
212 | iBinPair[0] = kMultBins; | |
213 | dBinsPair[0] = kMultBinLimits; | |
214 | axisTitlePair[0] = "kTPCITStracklet multiplicity"; | |
215 | } | |
216 | if(fEventClass == "Centrality"){ | |
683cbfb5 | 217 | // fine binning in case of vertex Z binning |
218 | if(fVertexBinning){ | |
219 | iBinSingle[0] = kNCentralityBinsVertex; | |
220 | dBinsSingle[0] = centralityBinsVertex; | |
221 | ||
222 | iBinPair[0] = kNCentralityBinsVertex; | |
223 | dBinsPair[0] = centralityBinsVertex; | |
224 | } | |
225 | else{ | |
226 | iBinSingle[0] = kNCentralityBins; | |
227 | dBinsSingle[0] = centralityBins; | |
228 | ||
229 | iBinPair[0] = kNCentralityBins; | |
230 | dBinsPair[0] = centralityBins; | |
231 | } | |
f0fb4ac1 | 232 | axisTitleSingle[0] = "Centrality percentile [%]"; |
f0fb4ac1 | 233 | axisTitlePair[0] = "Centrality percentile [%]"; |
234 | } | |
235 | if(fEventClass == "EventPlane"){ | |
236 | iBinSingle[0] = kNPsi2Bins; | |
237 | dBinsSingle[0] = psi2Bins; | |
238 | axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)"; | |
239 | iBinPair[0] = kNPsi2Bins; | |
240 | dBinsPair[0] = psi2Bins; | |
241 | axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)"; | |
242 | } | |
9afe3098 | 243 | |
683cbfb5 | 244 | // delta eta |
245 | const Int_t kNDeltaEtaBins = 80; | |
246 | const Int_t kNDeltaEtaBinsVertex = 40; | |
247 | Double_t deltaEtaBins[kNDeltaEtaBins+1]; | |
248 | Double_t deltaEtaBinsVertex[kNDeltaEtaBinsVertex+1]; | |
249 | for(Int_t i = 0; i < kNDeltaEtaBins+1; i++) | |
250 | deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins; | |
251 | for(Int_t i = 0; i < kNDeltaEtaBinsVertex+1; i++) | |
252 | deltaEtaBinsVertex[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBinsVertex; | |
253 | ||
254 | // coarse binning in case of vertex Z binning | |
255 | if(fVertexBinning){ | |
256 | iBinPair[1] = kNDeltaEtaBinsVertex; | |
257 | dBinsPair[1] = deltaEtaBinsVertex; | |
258 | } | |
259 | else{ | |
260 | iBinPair[1] = kNDeltaEtaBins; | |
261 | dBinsPair[1] = deltaEtaBins; | |
262 | } | |
263 | axisTitlePair[1] = "#Delta#eta"; | |
621329de | 264 | |
265 | // delta phi | |
9afe3098 | 266 | const Int_t kNDeltaPhiBins = 72; |
267 | Double_t deltaPhiBins[kNDeltaPhiBins+1]; | |
268 | for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){ | |
c8ad9635 | 269 | //deltaPhiBins[i] = -180.0 + i * 5.; |
270 | deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.; | |
9afe3098 | 271 | } |
621329de | 272 | iBinPair[2] = kNDeltaPhiBins; |
273 | dBinsPair[2] = deltaPhiBins; | |
c8ad9635 | 274 | axisTitlePair[2] = "#Delta#varphi (rad)"; |
9afe3098 | 275 | |
a38fd7f3 | 276 | // pt(trigger-associated) |
683cbfb5 | 277 | const Int_t kNPtBins = 16; |
a9f20288 | 278 | const Int_t kNPtBinsVertex = 6; |
683cbfb5 | 279 | Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.}; |
a9f20288 | 280 | Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0,15.0}; |
683cbfb5 | 281 | |
282 | // coarse binning in case of vertex Z binning | |
283 | if(fVertexBinning){ | |
284 | iBinSingle[1] = kNPtBinsVertex; | |
285 | dBinsSingle[1] = ptBinsVertex; | |
286 | ||
287 | iBinPair[3] = kNPtBinsVertex; | |
288 | dBinsPair[3] = ptBinsVertex; | |
289 | ||
290 | iBinPair[4] = kNPtBinsVertex; | |
291 | dBinsPair[4] = ptBinsVertex; | |
292 | } | |
293 | else{ | |
294 | iBinSingle[1] = kNPtBins; | |
295 | dBinsSingle[1] = ptBins; | |
296 | ||
297 | iBinPair[3] = kNPtBins; | |
298 | dBinsPair[3] = ptBins; | |
299 | ||
300 | iBinPair[4] = kNPtBins; | |
301 | dBinsPair[4] = ptBins; | |
302 | } | |
303 | ||
c8ad9635 | 304 | axisTitleSingle[1] = "p_{T,trig.} (GeV/c)"; |
683cbfb5 | 305 | axisTitlePair[3] = "p_{T,trig.} (GeV/c)"; |
306 | axisTitlePair[4] = "p_{T,assoc.} (GeV/c)"; | |
307 | ||
308 | // vertex Z | |
309 | const Int_t kNVertexZBins = 1; | |
310 | const Int_t kNVertexZBinsVertex = 9; | |
311 | Double_t vertexZBins[kNVertexZBins+1] = {-10., 10.}; | |
312 | Double_t vertexZBinsVertex[kNVertexZBinsVertex+1] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; | |
313 | ||
314 | // vertex Z binning or not | |
315 | if(fVertexBinning){ | |
316 | iBinSingle[2] = kNVertexZBinsVertex; | |
317 | dBinsSingle[2] = vertexZBinsVertex; | |
318 | ||
319 | iBinPair[5] = kNVertexZBinsVertex; | |
320 | dBinsPair[5] = vertexZBinsVertex; | |
321 | } | |
322 | else{ | |
323 | iBinSingle[2] = kNVertexZBins; | |
324 | dBinsSingle[2] = vertexZBins; | |
621329de | 325 | |
683cbfb5 | 326 | iBinPair[5] = kNVertexZBins; |
327 | dBinsPair[5] = vertexZBins; | |
328 | } | |
329 | ||
330 | axisTitleSingle[2] = "v_{Z} (cm)"; | |
331 | axisTitlePair[5] = "v_{Z} (cm)"; | |
a38fd7f3 | 332 | |
9afe3098 | 333 | |
334 | TString histName; | |
335 | //+ triggered particles | |
336 | histName = "fHistP"; | |
9fd4b54e | 337 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 338 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 339 | fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle); |
340 | for (Int_t j=0; j<kTrackVariablesSingle; j++) { | |
9afe3098 | 341 | fHistP->SetBinLimits(j, dBinsSingle[j]); |
342 | fHistP->SetVarTitle(j, axisTitleSingle[j]); | |
0879e280 | 343 | } |
9afe3098 | 344 | |
345 | //- triggered particles | |
346 | histName = "fHistN"; | |
9fd4b54e | 347 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 348 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 349 | fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle); |
350 | for (Int_t j=0; j<kTrackVariablesSingle; j++) { | |
9afe3098 | 351 | fHistN->SetBinLimits(j, dBinsSingle[j]); |
352 | fHistN->SetVarTitle(j, axisTitleSingle[j]); | |
0879e280 | 353 | } |
9afe3098 | 354 | |
355 | //+- pairs | |
356 | histName = "fHistPN"; | |
9fd4b54e | 357 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 358 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 359 | fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair); |
360 | for (Int_t j=0; j<kTrackVariablesPair; j++) { | |
9afe3098 | 361 | fHistPN->SetBinLimits(j, dBinsPair[j]); |
362 | fHistPN->SetVarTitle(j, axisTitlePair[j]); | |
0879e280 | 363 | } |
0879e280 | 364 | |
9afe3098 | 365 | //-+ pairs |
366 | histName = "fHistNP"; | |
9fd4b54e | 367 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 368 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 369 | fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair); |
370 | for (Int_t j=0; j<kTrackVariablesPair; j++) { | |
9afe3098 | 371 | fHistNP->SetBinLimits(j, dBinsPair[j]); |
372 | fHistNP->SetVarTitle(j, axisTitlePair[j]); | |
0879e280 | 373 | } |
0879e280 | 374 | |
9afe3098 | 375 | //++ pairs |
376 | histName = "fHistPP"; | |
9fd4b54e | 377 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 378 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 379 | fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair); |
380 | for (Int_t j=0; j<kTrackVariablesPair; j++) { | |
9afe3098 | 381 | fHistPP->SetBinLimits(j, dBinsPair[j]); |
382 | fHistPP->SetVarTitle(j, axisTitlePair[j]); | |
383 | } | |
384 | ||
385 | //-- pairs | |
386 | histName = "fHistNN"; | |
9fd4b54e | 387 | if(fShuffle) histName.Append("_shuffle"); |
9afe3098 | 388 | if(fCentralityId) histName += fCentralityId.Data(); |
9fd4b54e | 389 | fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair); |
390 | for (Int_t j=0; j<kTrackVariablesPair; j++) { | |
9afe3098 | 391 | fHistNN->SetBinLimits(j, dBinsPair[j]); |
392 | fHistNN->SetVarTitle(j, axisTitlePair[j]); | |
0879e280 | 393 | } |
f06d59b3 | 394 | AliInfo("Finished setting up the AliTHn"); |
73609a48 | 395 | |
396 | // QA histograms | |
c8ad9635 | 397 | fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi()); |
398 | fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi()); | |
03f5696b | 399 | fHistConversionbefore = new TH3D("fHistConversionbefore","before Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); |
400 | fHistConversionafter = new TH3D("fHistConversionafter","after Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
f0754617 | 401 | fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi()); |
402 | fHistResonancesBefore = new TH3D("fHistResonancesBefore","before resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
403 | fHistResonancesRho = new TH3D("fHistResonancesRho","after #rho resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
404 | fHistResonancesK0 = new TH3D("fHistResonancesK0","after #rho, K0 resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
405 | fHistResonancesLambda = new TH3D("fHistResonancesLambda","after #rho, K0, Lambda resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
41a6ae06 | 406 | fHistQbefore = new TH3D("fHistQbefore","before momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); |
407 | fHistQafter = new TH3D("fHistQafter","after momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5); | |
211b716d | 408 | |
409 | TH1::AddDirectory(oldStatus); | |
410 | ||
0879e280 | 411 | } |
412 | ||
413 | //____________________________________________________________________// | |
f06d59b3 | 414 | void AliBalancePsi::CalculateBalance(Double_t gReactionPlane, |
415 | TObjArray *particles, | |
73609a48 | 416 | TObjArray *particlesMixed, |
683cbfb5 | 417 | Float_t bSign, |
418 | Double_t kMultorCent, | |
419 | Double_t vertexZ) { | |
0879e280 | 420 | // Calculates the balance function |
421 | fAnalyzedEvents++; | |
f0fb4ac1 | 422 | |
0879e280 | 423 | // Initialize histograms if not done yet |
9afe3098 | 424 | if(!fHistPN){ |
0879e280 | 425 | AliWarning("Histograms not yet initialized! --> Will be done now"); |
426 | AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction"); | |
427 | InitHistograms(); | |
428 | } | |
429 | ||
9fd4b54e | 430 | Double_t trackVariablesSingle[kTrackVariablesSingle]; |
431 | Double_t trackVariablesPair[kTrackVariablesPair]; | |
f06d59b3 | 432 | |
433 | if (!particles){ | |
434 | AliWarning("particles TObjArray is NULL pointer --> return"); | |
435 | return; | |
436 | } | |
9afe3098 | 437 | |
f06d59b3 | 438 | // define end of particle loops |
439 | Int_t iMax = particles->GetEntriesFast(); | |
440 | Int_t jMax = iMax; | |
441 | if (particlesMixed) | |
442 | jMax = particlesMixed->GetEntriesFast(); | |
443 | ||
444 | // Eta() is extremely time consuming, therefore cache it for the inner loop here: | |
445 | TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles; | |
446 | ||
447 | TArrayF secondEta(jMax); | |
448 | TArrayF secondPhi(jMax); | |
449 | TArrayF secondPt(jMax); | |
c1847400 | 450 | TArrayS secondCharge(jMax); |
35aff0f3 | 451 | TArrayD secondCorrection(jMax); |
f06d59b3 | 452 | |
453 | for (Int_t i=0; i<jMax; i++){ | |
454 | secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta(); | |
455 | secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi(); | |
456 | secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt(); | |
c1847400 | 457 | secondCharge[i] = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge(); |
35aff0f3 | 458 | secondCorrection[i] = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction(); //==========================correction |
f06d59b3 | 459 | } |
460 | ||
2922bbe3 | 461 | //TLorenzVector implementation for resonances |
462 | TLorentzVector vectorMother, vectorDaughter[2]; | |
463 | TParticle pPion, pProton, pRho0, pK0s, pLambda; | |
464 | pPion.SetPdgCode(211); //pion | |
465 | pRho0.SetPdgCode(113); //rho0 | |
466 | pK0s.SetPdgCode(310); //K0s | |
467 | pProton.SetPdgCode(2212); //proton | |
468 | pLambda.SetPdgCode(3122); //Lambda | |
469 | Double_t gWidthForRho0 = 0.01; | |
470 | Double_t gWidthForK0s = 0.01; | |
471 | Double_t gWidthForLambda = 0.006; | |
472 | Double_t nSigmaRejection = 3.0; | |
473 | ||
f06d59b3 | 474 | // 1st particle loop |
73609a48 | 475 | for (Int_t i = 0; i < iMax; i++) { |
35aff0f3 | 476 | //AliVParticle* firstParticle = (AliVParticle*) particles->At(i); |
477 | AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction | |
6acdbcb2 | 478 | |
479 | // some optimization | |
480 | Float_t firstEta = firstParticle->Eta(); | |
481 | Float_t firstPhi = firstParticle->Phi(); | |
482 | Float_t firstPt = firstParticle->Pt(); | |
35aff0f3 | 483 | Float_t firstCorrection = firstParticle->Correction();//==========================correction |
484 | ||
6acdbcb2 | 485 | // Event plane (determine psi bin) |
486 | Double_t gPsiMinusPhi = 0.; | |
487 | Double_t gPsiMinusPhiBin = -10.; | |
488 | gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane); | |
489 | //in-plane | |
c8ad9635 | 490 | if((gPsiMinusPhi <= 7.5*TMath::DegToRad())|| |
491 | ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad()))) | |
6acdbcb2 | 492 | gPsiMinusPhiBin = 0.0; |
493 | //intermediate | |
c8ad9635 | 494 | else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))|| |
495 | ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))|| | |
496 | ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))|| | |
497 | ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad()))) | |
6acdbcb2 | 498 | gPsiMinusPhiBin = 1.0; |
499 | //out of plane | |
c8ad9635 | 500 | else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))|| |
501 | ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad()))) | |
6acdbcb2 | 502 | gPsiMinusPhiBin = 2.0; |
503 | //everything else | |
504 | else | |
505 | gPsiMinusPhiBin = 3.0; | |
506 | ||
c8ad9635 | 507 | fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi); |
73609a48 | 508 | |
509 | Short_t charge1 = (Short_t) firstParticle->Charge(); | |
6acdbcb2 | 510 | |
511 | trackVariablesSingle[0] = gPsiMinusPhiBin; | |
f0fb4ac1 | 512 | trackVariablesSingle[1] = firstPt; |
bb473a33 | 513 | if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent; |
683cbfb5 | 514 | trackVariablesSingle[2] = vertexZ; |
515 | ||
6acdbcb2 | 516 | |
517 | //fill single particle histograms | |
35aff0f3 | 518 | if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction |
519 | else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction | |
6acdbcb2 | 520 | |
c1847400 | 521 | // 2nd particle loop |
522 | for(Int_t j = 0; j < jMax; j++) { | |
523 | ||
524 | if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing) | |
525 | ||
d26c6047 | 526 | // pT,Assoc < pT,Trig |
527 | if(firstPt < secondPt[j]) continue; | |
528 | ||
c1847400 | 529 | Short_t charge2 = secondCharge[j]; |
9afe3098 | 530 | |
f0fb4ac1 | 531 | trackVariablesPair[0] = trackVariablesSingle[0]; |
6acdbcb2 | 532 | trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta |
533 | trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi | |
c8ad9635 | 534 | //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180 |
535 | //trackVariablesPair[2] -= 360.; | |
536 | //if (trackVariablesPair[2] < - 180.) | |
537 | //trackVariablesPair[2] += 360.; | |
538 | if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi | |
539 | trackVariablesPair[2] -= 2.*TMath::Pi(); | |
540 | if (trackVariablesPair[2] < - TMath::Pi()) | |
541 | trackVariablesPair[2] += 2.*TMath::Pi(); | |
542 | if (trackVariablesPair[2] < - TMath::Pi()/2.) | |
543 | trackVariablesPair[2] += 2.*TMath::Pi(); | |
9afe3098 | 544 | |
6acdbcb2 | 545 | trackVariablesPair[3] = firstPt; // pt trigger |
546 | trackVariablesPair[4] = secondPt[j]; // pt | |
683cbfb5 | 547 | trackVariablesPair[5] = vertexZ; // z of the primary vertex |
2922bbe3 | 548 | |
549 | //Exclude resonances for the calculation of pairs by looking | |
550 | //at the invariant mass and not considering the pairs that | |
551 | //fall within 3sigma from the mass peak of: rho0, K0s, Lambda | |
552 | if(fResonancesCut) { | |
749de081 | 553 | if (charge1 * charge2 < 0) { |
554 | ||
555 | //rho0 | |
556 | vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass()); | |
557 | vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass()); | |
558 | vectorMother = vectorDaughter[0] + vectorDaughter[1]; | |
559 | fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M()); | |
560 | if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0) | |
561 | continue; | |
562 | fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M()); | |
563 | ||
564 | //K0s | |
565 | if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s) | |
566 | continue; | |
567 | fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M()); | |
568 | ||
569 | ||
570 | //Lambda | |
571 | vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass()); | |
572 | vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass()); | |
573 | vectorMother = vectorDaughter[0] + vectorDaughter[1]; | |
574 | if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda) | |
575 | continue; | |
576 | ||
577 | vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass()); | |
578 | vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass()); | |
579 | vectorMother = vectorDaughter[0] + vectorDaughter[1]; | |
580 | if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda) | |
581 | continue; | |
582 | fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M()); | |
583 | ||
584 | }//unlike-sign only | |
2922bbe3 | 585 | }//resonance cut |
73609a48 | 586 | |
587 | // HBT like cut | |
da8b2d30 | 588 | if(fHBTCut){ // VERSION 3 (all pairs) |
589 | //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS) | |
73609a48 | 590 | //if( dphi < 3 || deta < 0.01 ){ // VERSION 1 |
591 | // continue; | |
592 | ||
593 | Double_t deta = firstEta - secondEta[j]; | |
594 | Double_t dphi = firstPhi - secondPhi[j]; | |
595 | // VERSION 2 (Taken from DPhiCorrelations) | |
596 | // the variables & cuthave been developed by the HBT group | |
597 | // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700 | |
598 | fHistHBTbefore->Fill(deta,dphi); | |
599 | ||
600 | // optimization | |
d9eb282c | 601 | if (TMath::Abs(deta) < fHBTCutValue * 2.5 * 3) //fHBTCutValue = 0.02 [default for dphicorrelations] |
73609a48 | 602 | { |
603 | // phi in rad | |
c8ad9635 | 604 | //Float_t phi1rad = firstPhi*TMath::DegToRad(); |
605 | //Float_t phi2rad = secondPhi[j]*TMath::DegToRad(); | |
606 | Float_t phi1rad = firstPhi; | |
607 | Float_t phi2rad = secondPhi[j]; | |
73609a48 | 608 | |
609 | // check first boundaries to see if is worth to loop and find the minimum | |
610 | Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign); | |
611 | Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign); | |
612 | ||
613 | const Float_t kLimit = 0.02 * 3; | |
614 | ||
615 | Float_t dphistarminabs = 1e5; | |
616 | Float_t dphistarmin = 1e5; | |
617 | ||
618 | if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) { | |
619 | for (Double_t rad=0.8; rad<2.51; rad+=0.01) { | |
620 | Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign); | |
621 | Float_t dphistarabs = TMath::Abs(dphistar); | |
622 | ||
623 | if (dphistarabs < dphistarminabs) { | |
624 | dphistarmin = dphistar; | |
625 | dphistarminabs = dphistarabs; | |
626 | } | |
627 | } | |
628 | ||
629 | if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) { | |
630 | //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign)); | |
631 | continue; | |
632 | } | |
633 | } | |
634 | } | |
635 | fHistHBTafter->Fill(deta,dphi); | |
636 | }//HBT cut | |
637 | ||
638 | // conversions | |
639 | if(fConversionCut) { | |
640 | if (charge1 * charge2 < 0) { | |
641 | Double_t deta = firstEta - secondEta[j]; | |
642 | Double_t dphi = firstPhi - secondPhi[j]; | |
73609a48 | 643 | |
644 | Float_t m0 = 0.510e-3; | |
645 | Float_t tantheta1 = 1e10; | |
646 | ||
647 | // phi in rad | |
c8ad9635 | 648 | //Float_t phi1rad = firstPhi*TMath::DegToRad(); |
649 | //Float_t phi2rad = secondPhi[j]*TMath::DegToRad(); | |
650 | Float_t phi1rad = firstPhi; | |
651 | Float_t phi2rad = secondPhi[j]; | |
73609a48 | 652 | |
653 | if (firstEta < -1e-10 || firstEta > 1e-10) | |
654 | tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta)); | |
655 | ||
656 | Float_t tantheta2 = 1e10; | |
657 | if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10) | |
658 | tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j])); | |
659 | ||
660 | Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1); | |
661 | Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2); | |
662 | ||
663 | Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) ); | |
03f5696b | 664 | |
665 | fHistConversionbefore->Fill(deta,dphi,masssqu); | |
73609a48 | 666 | |
a9737921 | 667 | if (masssqu < fInvMassCutConversion*fInvMassCutConversion){ |
73609a48 | 668 | //AliInfo(Form("Conversion: Removed track pair %d %d with [[%f %f] %f %f] %d %d <- %f %f %f %f %f %f ", i, j, deta, dphi, masssqu, charge1, charge2,eta1,eta2,phi1,phi2,pt1,pt2)); |
669 | continue; | |
670 | } | |
03f5696b | 671 | fHistConversionafter->Fill(deta,dphi,masssqu); |
73609a48 | 672 | } |
673 | }//conversion cut | |
5a861dc6 | 674 | |
675 | // momentum difference cut - suppress femtoscopic effects | |
676 | if(fQCut){ | |
677 | ||
6fa567bd | 678 | //Double_t ptMin = 0.1; //const for the time being (should be changeable later on) |
41a6ae06 | 679 | Double_t ptDifference = TMath::Abs( firstPt - secondPt[j]); |
5a861dc6 | 680 | |
41a6ae06 | 681 | fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference); |
6fa567bd | 682 | if(ptDifference < fDeltaPtMin) continue; |
41a6ae06 | 683 | fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference); |
5a861dc6 | 684 | |
685 | } | |
686 | ||
35aff0f3 | 687 | if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction |
688 | else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction | |
689 | else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction | |
690 | else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction | |
8795e993 | 691 | else { |
692 | //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2)); | |
693 | continue; | |
694 | } | |
6acdbcb2 | 695 | }//end of 2nd particle loop |
696 | }//end of 1st particle loop | |
0879e280 | 697 | } |
698 | ||
0879e280 | 699 | //____________________________________________________________________// |
9afe3098 | 700 | TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle, |
701 | Int_t iVariablePair, | |
9afe3098 | 702 | Double_t psiMin, |
6acdbcb2 | 703 | Double_t psiMax, |
20006629 | 704 | Double_t vertexZMin, |
705 | Double_t vertexZMax, | |
6acdbcb2 | 706 | Double_t ptTriggerMin, |
707 | Double_t ptTriggerMax, | |
708 | Double_t ptAssociatedMin, | |
709 | Double_t ptAssociatedMax) { | |
9afe3098 | 710 | //Returns the BF histogram, extracted from the 6 AliTHn objects |
0879e280 | 711 | //(private members) of the AliBalancePsi class. |
621329de | 712 | //iVariableSingle: 0(phi-Psi), 1(pt-trigger) |
713 | //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated | |
622473b9 | 714 | |
715 | // security checks | |
b0d53e32 | 716 | if(psiMin > psiMax-0.00001){ |
717 | AliError("psiMax <= psiMin"); | |
718 | return NULL; | |
719 | } | |
20006629 | 720 | if(vertexZMin > vertexZMax-0.00001){ |
721 | AliError("vertexZMax <= vertexZMin"); | |
722 | return NULL; | |
723 | } | |
b0d53e32 | 724 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
725 | AliError("ptTriggerMax <= ptTriggerMin"); | |
726 | return NULL; | |
727 | } | |
728 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
729 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
730 | return NULL; | |
731 | } | |
622473b9 | 732 | |
9afe3098 | 733 | // Psi_2 |
622473b9 | 734 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
735 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
736 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
737 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
738 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
739 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
9afe3098 | 740 | |
20006629 | 741 | // Vz |
742 | if(fVertexBinning){ | |
743 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
744 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
745 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
746 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
747 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
748 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
749 | } | |
750 | ||
6acdbcb2 | 751 | // pt trigger |
752 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 753 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
754 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
755 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
756 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
757 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
758 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
6acdbcb2 | 759 | } |
760 | ||
761 | // pt associated | |
762 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { | |
622473b9 | 763 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
764 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
765 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
766 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
6acdbcb2 | 767 | } |
768 | ||
a38fd7f3 | 769 | //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)); |
770 | ||
9afe3098 | 771 | // Project into the wanted space (1st: analysis step, 2nd: axis) |
20006629 | 772 | TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); // |
773 | TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); // | |
774 | TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); // | |
775 | TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); // | |
776 | TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); // | |
777 | TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); // | |
9afe3098 | 778 | |
779 | TH1D *gHistBalanceFunctionHistogram = 0x0; | |
780 | if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) { | |
781 | gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone(); | |
782 | gHistBalanceFunctionHistogram->Reset(); | |
783 | ||
784 | switch(iVariablePair) { | |
621329de | 785 | case 1: |
c8ad9635 | 786 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta"); |
787 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)"); | |
9afe3098 | 788 | break; |
621329de | 789 | case 2: |
c8ad9635 | 790 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)"); |
791 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)"); | |
9afe3098 | 792 | break; |
9afe3098 | 793 | default: |
794 | break; | |
795 | } | |
0879e280 | 796 | |
0879e280 | 797 | hTemp1->Sumw2(); |
798 | hTemp2->Sumw2(); | |
799 | hTemp3->Sumw2(); | |
800 | hTemp4->Sumw2(); | |
a38fd7f3 | 801 | hTemp1->Add(hTemp3,-1.); |
0879e280 | 802 | hTemp1->Scale(1./hTemp5->GetEntries()); |
a38fd7f3 | 803 | hTemp2->Add(hTemp4,-1.); |
0879e280 | 804 | hTemp2->Scale(1./hTemp6->GetEntries()); |
805 | gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.); | |
9afe3098 | 806 | gHistBalanceFunctionHistogram->Scale(0.5); |
faa03677 | 807 | |
808 | //normalize to bin width | |
a9f20288 | 809 | gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1))); |
0879e280 | 810 | } |
811 | ||
0879e280 | 812 | return gHistBalanceFunctionHistogram; |
813 | } | |
a38fd7f3 | 814 | |
f0c5040c | 815 | //____________________________________________________________________// |
816 | TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle, | |
817 | Int_t iVariablePair, | |
818 | Double_t psiMin, | |
819 | Double_t psiMax, | |
20006629 | 820 | Double_t vertexZMin, |
821 | Double_t vertexZMax, | |
f0c5040c | 822 | Double_t ptTriggerMin, |
823 | Double_t ptTriggerMax, | |
824 | Double_t ptAssociatedMin, | |
825 | Double_t ptAssociatedMax, | |
826 | AliBalancePsi *bfMix) { | |
827 | //Returns the BF histogram, extracted from the 6 AliTHn objects | |
828 | //after dividing each correlation function by the Event Mixing one | |
829 | //(private members) of the AliBalancePsi class. | |
830 | //iVariableSingle: 0(phi-Psi), 1(pt-trigger) | |
831 | //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated | |
622473b9 | 832 | |
833 | // security checks | |
b0d53e32 | 834 | if(psiMin > psiMax-0.00001){ |
835 | AliError("psiMax <= psiMin"); | |
836 | return NULL; | |
837 | } | |
20006629 | 838 | if(vertexZMin > vertexZMax-0.00001){ |
839 | AliError("vertexZMax <= vertexZMin"); | |
840 | return NULL; | |
841 | } | |
b0d53e32 | 842 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
843 | AliError("ptTriggerMax <= ptTriggerMin"); | |
844 | return NULL; | |
845 | } | |
846 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
847 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
848 | return NULL; | |
849 | } | |
20006629 | 850 | |
34616eda | 851 | // pt trigger (fixed for all vertices/psi/centralities) |
f0c5040c | 852 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 853 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
854 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
855 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
856 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
857 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
858 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f0c5040c | 859 | } |
860 | ||
34616eda | 861 | // pt associated (fixed for all vertices/psi/centralities) |
f0c5040c | 862 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 863 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
864 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
865 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
866 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f0c5040c | 867 | } |
868 | ||
869 | // ============================================================================================ | |
870 | // the same for event mixing | |
871 | AliTHn *fHistPMix = bfMix->GetHistNp(); | |
872 | AliTHn *fHistNMix = bfMix->GetHistNn(); | |
873 | AliTHn *fHistPNMix = bfMix->GetHistNpn(); | |
874 | AliTHn *fHistNPMix = bfMix->GetHistNnp(); | |
875 | AliTHn *fHistPPMix = bfMix->GetHistNpp(); | |
876 | AliTHn *fHistNNMix = bfMix->GetHistNnn(); | |
877 | ||
34616eda | 878 | // pt trigger (fixed for all vertices/psi/centralities) |
f0c5040c | 879 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 880 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
881 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
882 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
883 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
884 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
885 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f0c5040c | 886 | } |
887 | ||
34616eda | 888 | // pt associated (fixed for all vertices/psi/centralities) |
f0c5040c | 889 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 890 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
891 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
892 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
893 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f0c5040c | 894 | } |
34616eda | 895 | |
f0c5040c | 896 | // ============================================================================================ |
34616eda | 897 | // ranges (in bins) for vertexZ and centrality (psi) |
898 | Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin); | |
899 | Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001); | |
a4e2c68f | 900 | Int_t binVertexMin = 0; |
901 | Int_t binVertexMax = 0; | |
902 | if(fVertexBinning){ | |
903 | binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin); | |
904 | binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001); | |
905 | } | |
34616eda | 906 | Double_t binPsiLowEdge = 0.; |
a4e2c68f | 907 | Double_t binPsiUpEdge = 1000.; |
34616eda | 908 | Double_t binVertexLowEdge = 0.; |
a4e2c68f | 909 | Double_t binVertexUpEdge = 1000.; |
34616eda | 910 | |
911 | TH1D* h1 = NULL; | |
912 | TH1D* h2 = NULL; | |
913 | TH1D* h3 = NULL; | |
914 | TH1D* h4 = NULL; | |
915 | ||
916 | // loop over all bins | |
917 | for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){ | |
918 | for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){ | |
919 | ||
920 | cout<<"In the balance function (1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl; | |
f0c5040c | 921 | |
34616eda | 922 | // determine the bin edges for this bin |
923 | binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi); | |
924 | binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi); | |
a4e2c68f | 925 | if(fVertexBinning){ |
926 | binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex); | |
927 | binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex); | |
928 | } | |
34616eda | 929 | |
930 | // Psi_2 | |
931 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
932 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
933 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
934 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
935 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
936 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
937 | ||
938 | // Vz | |
939 | if(fVertexBinning){ | |
940 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
941 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
942 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
943 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
944 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
945 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
946 | } | |
f0c5040c | 947 | |
34616eda | 948 | // ============================================================================================ |
949 | // the same for event mixing | |
950 | ||
951 | // Psi_2 | |
952 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
953 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
954 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
955 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
956 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
957 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
958 | ||
959 | // Vz | |
960 | if(fVertexBinning){ | |
961 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
962 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
963 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
964 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
965 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
966 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
967 | } | |
968 | ||
969 | // ============================================================================================ | |
f0c5040c | 970 | |
34616eda | 971 | //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)); |
972 | ||
973 | // Project into the wanted space (1st: analysis step, 2nd: axis) | |
974 | TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); | |
975 | TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); | |
976 | TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); | |
977 | TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); | |
978 | TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); | |
979 | TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); | |
980 | ||
981 | // ============================================================================================ | |
982 | // the same for event mixing | |
983 | TH1D* hTemp1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair); | |
984 | TH1D* hTemp2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair); | |
985 | TH1D* hTemp3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair); | |
986 | TH1D* hTemp4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair); | |
987 | TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle); | |
988 | TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle); | |
989 | // ============================================================================================ | |
990 | ||
991 | hTemp1->Sumw2(); | |
992 | hTemp2->Sumw2(); | |
993 | hTemp3->Sumw2(); | |
994 | hTemp4->Sumw2(); | |
995 | hTemp1Mix->Sumw2(); | |
996 | hTemp2Mix->Sumw2(); | |
997 | hTemp3Mix->Sumw2(); | |
998 | hTemp4Mix->Sumw2(); | |
999 | ||
1000 | hTemp1->Scale(1./hTemp5->GetEntries()); | |
1001 | hTemp3->Scale(1./hTemp5->GetEntries()); | |
1002 | hTemp2->Scale(1./hTemp6->GetEntries()); | |
1003 | hTemp4->Scale(1./hTemp6->GetEntries()); | |
271ceae3 | 1004 | |
1005 | // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method | |
1006 | // does not work here, so normalize also to trigger particles | |
1007 | // --> careful: gives different integrals then as with full 2D method | |
34616eda | 1008 | hTemp1Mix->Scale(1./hTemp5Mix->GetEntries()); |
1009 | hTemp3Mix->Scale(1./hTemp5Mix->GetEntries()); | |
1010 | hTemp2Mix->Scale(1./hTemp6Mix->GetEntries()); | |
1011 | hTemp4Mix->Scale(1./hTemp6Mix->GetEntries()); | |
271ceae3 | 1012 | |
34616eda | 1013 | hTemp1->Divide(hTemp1Mix); |
1014 | hTemp2->Divide(hTemp2Mix); | |
1015 | hTemp3->Divide(hTemp3Mix); | |
1016 | hTemp4->Divide(hTemp4Mix); | |
1017 | ||
1018 | // for the first: clone | |
1019 | if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){ | |
1020 | h1 = (TH1D*)hTemp1->Clone(); | |
1021 | h2 = (TH1D*)hTemp2->Clone(); | |
1022 | h3 = (TH1D*)hTemp3->Clone(); | |
1023 | h4 = (TH1D*)hTemp4->Clone(); | |
1024 | } | |
1025 | else{ // otherwise: add for averaging | |
1026 | h1->Add(hTemp1); | |
1027 | h2->Add(hTemp2); | |
1028 | h3->Add(hTemp3); | |
1029 | h4->Add(hTemp4); | |
1030 | } | |
1031 | ||
1032 | } | |
1033 | } | |
f0c5040c | 1034 | |
1035 | TH1D *gHistBalanceFunctionHistogram = 0x0; | |
34616eda | 1036 | if((h1)&&(h2)&&(h3)&&(h4)) { |
1037 | ||
1038 | // average over number of bins nbinsVertex * nbinsPsi | |
1039 | h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1040 | h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1041 | h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1042 | h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1043 | ||
1044 | gHistBalanceFunctionHistogram = (TH1D*)h1->Clone(); | |
f0c5040c | 1045 | gHistBalanceFunctionHistogram->Reset(); |
1046 | ||
1047 | switch(iVariablePair) { | |
1048 | case 1: | |
1049 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta"); | |
1050 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)"); | |
1051 | break; | |
1052 | case 2: | |
1053 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)"); | |
1054 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)"); | |
1055 | break; | |
1056 | default: | |
1057 | break; | |
1058 | } | |
1059 | ||
34616eda | 1060 | h1->Add(h3,-1.); |
1061 | h2->Add(h4,-1.); | |
f0c5040c | 1062 | |
34616eda | 1063 | gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.); |
f0c5040c | 1064 | gHistBalanceFunctionHistogram->Scale(0.5); |
faa03677 | 1065 | |
1066 | //normalize to bin width | |
a9f20288 | 1067 | gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1))); |
f0c5040c | 1068 | } |
1069 | ||
1070 | return gHistBalanceFunctionHistogram; | |
1071 | } | |
1072 | ||
a38fd7f3 | 1073 | //____________________________________________________________________// |
621329de | 1074 | TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin, |
6acdbcb2 | 1075 | Double_t psiMax, |
20006629 | 1076 | Double_t vertexZMin, |
1077 | Double_t vertexZMax, | |
6acdbcb2 | 1078 | Double_t ptTriggerMin, |
1079 | Double_t ptTriggerMax, | |
1080 | Double_t ptAssociatedMin, | |
1081 | Double_t ptAssociatedMax) { | |
621329de | 1082 | //Returns the BF histogram in Delta eta vs Delta phi, |
1083 | //extracted from the 6 AliTHn objects | |
1084 | //(private members) of the AliBalancePsi class. | |
1085 | //iVariableSingle: 0(phi-Psi), 1(pt-trigger) | |
1086 | //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated | |
622473b9 | 1087 | |
1088 | // security checks | |
b0d53e32 | 1089 | if(psiMin > psiMax-0.00001){ |
1090 | AliError("psiMax <= psiMin"); | |
1091 | return NULL; | |
1092 | } | |
20006629 | 1093 | if(vertexZMin > vertexZMax-0.00001){ |
1094 | AliError("vertexZMax <= vertexZMin"); | |
1095 | return NULL; | |
1096 | } | |
b0d53e32 | 1097 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
1098 | AliError("ptTriggerMax <= ptTriggerMin"); | |
1099 | return NULL; | |
1100 | } | |
1101 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
1102 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
1103 | return NULL; | |
1104 | } | |
622473b9 | 1105 | |
621329de | 1106 | TString histName = "gHistBalanceFunctionHistogram2D"; |
1107 | ||
1108 | // Psi_2 | |
622473b9 | 1109 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
1110 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1111 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1112 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1113 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1114 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
621329de | 1115 | |
20006629 | 1116 | // Vz |
1117 | if(fVertexBinning){ | |
1118 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1119 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1120 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1121 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1122 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1123 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1124 | } | |
1125 | ||
6acdbcb2 | 1126 | // pt trigger |
1127 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 1128 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1129 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1130 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1131 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1132 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1133 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
6acdbcb2 | 1134 | } |
1135 | ||
1136 | // pt associated | |
1137 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { | |
622473b9 | 1138 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
1139 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1140 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1141 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
6acdbcb2 | 1142 | } |
1143 | ||
1144 | //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0))); | |
621329de | 1145 | |
1146 | // Project into the wanted space (1st: analysis step, 2nd: axis) | |
1147 | TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2); | |
1148 | TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2); | |
1149 | TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2); | |
1150 | TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2); | |
1151 | TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1); | |
1152 | TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1); | |
1153 | ||
1154 | TH2D *gHistBalanceFunctionHistogram = 0x0; | |
1155 | if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) { | |
1156 | gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone(); | |
1157 | gHistBalanceFunctionHistogram->Reset(); | |
c8ad9635 | 1158 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta"); |
1159 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)"); | |
1160 | gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)"); | |
621329de | 1161 | |
1162 | hTemp1->Sumw2(); | |
1163 | hTemp2->Sumw2(); | |
1164 | hTemp3->Sumw2(); | |
1165 | hTemp4->Sumw2(); | |
1166 | hTemp1->Add(hTemp3,-1.); | |
1167 | hTemp1->Scale(1./hTemp5->GetEntries()); | |
1168 | hTemp2->Add(hTemp4,-1.); | |
1169 | hTemp2->Scale(1./hTemp6->GetEntries()); | |
1170 | gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.); | |
1171 | gHistBalanceFunctionHistogram->Scale(0.5); | |
faa03677 | 1172 | |
1173 | //normalize to bin width | |
1174 | gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1))); | |
621329de | 1175 | } |
1176 | ||
1177 | return gHistBalanceFunctionHistogram; | |
1178 | } | |
1179 | ||
f0c5040c | 1180 | //____________________________________________________________________// |
1181 | TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin, | |
1182 | Double_t psiMax, | |
20006629 | 1183 | Double_t vertexZMin, |
1184 | Double_t vertexZMax, | |
f0c5040c | 1185 | Double_t ptTriggerMin, |
1186 | Double_t ptTriggerMax, | |
1187 | Double_t ptAssociatedMin, | |
1188 | Double_t ptAssociatedMax, | |
1189 | AliBalancePsi *bfMix) { | |
1190 | //Returns the BF histogram in Delta eta vs Delta phi, | |
1191 | //after dividing each correlation function by the Event Mixing one | |
1192 | //extracted from the 6 AliTHn objects | |
1193 | //(private members) of the AliBalancePsi class. | |
1194 | //iVariableSingle: 0(phi-Psi), 1(pt-trigger) | |
1195 | //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated | |
1196 | TString histName = "gHistBalanceFunctionHistogram2D"; | |
1197 | ||
1198 | if(!bfMix){ | |
1199 | AliError("balance function object for event mixing not available"); | |
1200 | return NULL; | |
1201 | } | |
1202 | ||
622473b9 | 1203 | // security checks |
b0d53e32 | 1204 | if(psiMin > psiMax-0.00001){ |
1205 | AliError("psiMax <= psiMin"); | |
1206 | return NULL; | |
1207 | } | |
20006629 | 1208 | if(vertexZMin > vertexZMax-0.00001){ |
1209 | AliError("vertexZMax <= vertexZMin"); | |
1210 | return NULL; | |
1211 | } | |
b0d53e32 | 1212 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
1213 | AliError("ptTriggerMax <= ptTriggerMin"); | |
1214 | return NULL; | |
1215 | } | |
1216 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
1217 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
1218 | return NULL; | |
1219 | } | |
622473b9 | 1220 | |
34616eda | 1221 | // pt trigger (fixed for all vertices/psi/centralities) |
f0c5040c | 1222 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 1223 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1224 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1225 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1226 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1227 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1228 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f0c5040c | 1229 | } |
1230 | ||
34616eda | 1231 | // pt associated (fixed for all vertices/psi/centralities) |
f0c5040c | 1232 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 1233 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
1234 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1235 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1236 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f0c5040c | 1237 | } |
1238 | ||
f0c5040c | 1239 | // ============================================================================================ |
1240 | // the same for event mixing | |
1241 | AliTHn *fHistPMix = bfMix->GetHistNp(); | |
1242 | AliTHn *fHistNMix = bfMix->GetHistNn(); | |
1243 | AliTHn *fHistPNMix = bfMix->GetHistNpn(); | |
1244 | AliTHn *fHistNPMix = bfMix->GetHistNnp(); | |
1245 | AliTHn *fHistPPMix = bfMix->GetHistNpp(); | |
1246 | AliTHn *fHistNNMix = bfMix->GetHistNnn(); | |
1247 | ||
34616eda | 1248 | // pt trigger (fixed for all vertices/psi/centralities) |
f0c5040c | 1249 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 1250 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1251 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1252 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1253 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1254 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1255 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f0c5040c | 1256 | } |
1257 | ||
34616eda | 1258 | // pt associated (fixed for all vertices/psi/centralities) |
f0c5040c | 1259 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 1260 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
1261 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1262 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1263 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f0c5040c | 1264 | } |
34616eda | 1265 | |
f0c5040c | 1266 | // ============================================================================================ |
34616eda | 1267 | // ranges (in bins) for vertexZ and centrality (psi) |
1268 | Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin); | |
1269 | Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001); | |
a4e2c68f | 1270 | Int_t binVertexMin = 0; |
1271 | Int_t binVertexMax = 0; | |
1272 | if(fVertexBinning){ | |
1273 | binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin); | |
1274 | binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001); | |
1275 | } | |
34616eda | 1276 | Double_t binPsiLowEdge = 0.; |
1277 | Double_t binPsiUpEdge = 0.; | |
1278 | Double_t binVertexLowEdge = 0.; | |
1279 | Double_t binVertexUpEdge = 0.; | |
1280 | ||
1281 | TH2D* h1 = NULL; | |
1282 | TH2D* h2 = NULL; | |
1283 | TH2D* h3 = NULL; | |
1284 | TH2D* h4 = NULL; | |
1285 | ||
1286 | // loop over all bins | |
1287 | for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){ | |
1288 | for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){ | |
1289 | ||
1290 | cout<<"In the balance function (2D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl; | |
1291 | ||
1292 | // determine the bin edges for this bin | |
1293 | binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi); | |
1294 | binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi); | |
a4e2c68f | 1295 | if(fVertexBinning){ |
1296 | binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex); | |
1297 | binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex); | |
1298 | } | |
34616eda | 1299 | |
1300 | // Psi_2 | |
1301 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1302 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1303 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1304 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1305 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1306 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1307 | ||
1308 | // Vz | |
1309 | if(fVertexBinning){ | |
1310 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1311 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1312 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1313 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1314 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1315 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1316 | } | |
f0c5040c | 1317 | |
1318 | ||
f0c5040c | 1319 | |
34616eda | 1320 | // ============================================================================================ |
1321 | // the same for event mixing | |
1322 | ||
1323 | // Psi_2 | |
1324 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1325 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1326 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1327 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1328 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1329 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001); | |
1330 | ||
1331 | // Vz | |
1332 | if(fVertexBinning){ | |
1333 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1334 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1335 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1336 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1337 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1338 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001); | |
1339 | } | |
1340 | ||
1341 | // ============================================================================================ | |
1342 | ||
1343 | ||
1344 | //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0))); | |
1345 | ||
1346 | // Project into the wanted space (1st: analysis step, 2nd: axis) | |
1347 | TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2); | |
1348 | TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2); | |
1349 | TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2); | |
1350 | TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2); | |
1351 | TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1); | |
1352 | TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1); | |
1353 | ||
1354 | // ============================================================================================ | |
1355 | // the same for event mixing | |
1356 | TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2); | |
1357 | TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2); | |
1358 | TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2); | |
1359 | TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2); | |
d631c24b | 1360 | // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1); |
1361 | // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1); | |
34616eda | 1362 | // ============================================================================================ |
1363 | ||
1364 | hTemp1->Sumw2(); | |
1365 | hTemp2->Sumw2(); | |
1366 | hTemp3->Sumw2(); | |
1367 | hTemp4->Sumw2(); | |
1368 | hTemp1Mix->Sumw2(); | |
1369 | hTemp2Mix->Sumw2(); | |
1370 | hTemp3Mix->Sumw2(); | |
1371 | hTemp4Mix->Sumw2(); | |
1372 | ||
1373 | hTemp1->Scale(1./hTemp5->GetEntries()); | |
1374 | hTemp3->Scale(1./hTemp5->GetEntries()); | |
1375 | hTemp2->Scale(1./hTemp6->GetEntries()); | |
1376 | hTemp4->Scale(1./hTemp6->GetEntries()); | |
271ceae3 | 1377 | |
1378 | // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method | |
1379 | Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX()); | |
1380 | mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1381 | hTemp1Mix->Scale(1./mixedNorm1); | |
1382 | Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX()); | |
1383 | mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1384 | hTemp2Mix->Scale(1./mixedNorm2); | |
1385 | Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX()); | |
1386 | mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1387 | hTemp3Mix->Scale(1./mixedNorm3); | |
1388 | Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX()); | |
1389 | mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1390 | hTemp4Mix->Scale(1./mixedNorm4); | |
34616eda | 1391 | |
1392 | hTemp1->Divide(hTemp1Mix); | |
1393 | hTemp2->Divide(hTemp2Mix); | |
1394 | hTemp3->Divide(hTemp3Mix); | |
1395 | hTemp4->Divide(hTemp4Mix); | |
1396 | ||
1397 | // for the first: clone | |
1398 | if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){ | |
1399 | h1 = (TH2D*)hTemp1->Clone(); | |
1400 | h2 = (TH2D*)hTemp2->Clone(); | |
1401 | h3 = (TH2D*)hTemp3->Clone(); | |
1402 | h4 = (TH2D*)hTemp4->Clone(); | |
1403 | } | |
1404 | else{ // otherwise: add for averaging | |
1405 | h1->Add(hTemp1); | |
1406 | h2->Add(hTemp2); | |
1407 | h3->Add(hTemp3); | |
1408 | h4->Add(hTemp4); | |
1409 | } | |
1410 | } | |
1411 | } | |
1412 | ||
1413 | TH2D *gHistBalanceFunctionHistogram = 0x0; | |
1414 | if((h1)&&(h2)&&(h3)&&(h4)) { | |
f0c5040c | 1415 | |
34616eda | 1416 | // average over number of bins nbinsVertex * nbinsPsi |
1417 | h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1418 | h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1419 | h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1420 | h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
f0c5040c | 1421 | |
34616eda | 1422 | gHistBalanceFunctionHistogram = (TH2D*)h1->Clone(); |
f0c5040c | 1423 | gHistBalanceFunctionHistogram->Reset(); |
1424 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta"); | |
1425 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)"); | |
1426 | gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)"); | |
f0c5040c | 1427 | |
34616eda | 1428 | h1->Add(h3,-1.); |
1429 | h2->Add(h4,-1.); | |
1430 | ||
1431 | gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.); | |
f0c5040c | 1432 | gHistBalanceFunctionHistogram->Scale(0.5); |
34616eda | 1433 | |
faa03677 | 1434 | //normalize to bin width |
1435 | gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1))); | |
f0c5040c | 1436 | } |
34616eda | 1437 | |
f0c5040c | 1438 | return gHistBalanceFunctionHistogram; |
1439 | } | |
1440 | ||
f2e8af26 | 1441 | //____________________________________________________________________// |
1442 | TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi, | |
1443 | Double_t psiMin, | |
1444 | Double_t psiMax, | |
20006629 | 1445 | Double_t vertexZMin, |
1446 | Double_t vertexZMax, | |
f2e8af26 | 1447 | Double_t ptTriggerMin, |
1448 | Double_t ptTriggerMax, | |
1449 | Double_t ptAssociatedMin, | |
1450 | Double_t ptAssociatedMax, | |
1451 | AliBalancePsi *bfMix) { | |
1452 | //Returns the BF histogram in Delta eta OR Delta phi, | |
1453 | //after dividing each correlation function by the Event Mixing one | |
1454 | // (But the division is done here in 2D, this was basically done to check the Integral) | |
1455 | //extracted from the 6 AliTHn objects | |
1456 | //(private members) of the AliBalancePsi class. | |
1457 | //iVariableSingle: 0(phi-Psi), 1(pt-trigger) | |
1458 | //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated | |
1459 | TString histName = "gHistBalanceFunctionHistogram1D"; | |
1460 | ||
1461 | if(!bfMix){ | |
1462 | AliError("balance function object for event mixing not available"); | |
1463 | return NULL; | |
1464 | } | |
1465 | ||
622473b9 | 1466 | // security checks |
b0d53e32 | 1467 | if(psiMin > psiMax-0.00001){ |
1468 | AliError("psiMax <= psiMin"); | |
1469 | return NULL; | |
1470 | } | |
20006629 | 1471 | if(vertexZMin > vertexZMax-0.00001){ |
1472 | AliError("vertexZMax <= vertexZMin"); | |
1473 | return NULL; | |
1474 | } | |
b0d53e32 | 1475 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
1476 | AliError("ptTriggerMax <= ptTriggerMin"); | |
1477 | return NULL; | |
1478 | } | |
1479 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
1480 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
1481 | return NULL; | |
1482 | } | |
622473b9 | 1483 | |
34616eda | 1484 | // pt trigger (fixed for all vertices/psi/centralities) |
f2e8af26 | 1485 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 1486 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1487 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1488 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1489 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1490 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1491 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f2e8af26 | 1492 | } |
1493 | ||
34616eda | 1494 | // pt associated (fixed for all vertices/psi/centralities) |
f2e8af26 | 1495 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 1496 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
1497 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1498 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1499 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f2e8af26 | 1500 | } |
1501 | ||
f2e8af26 | 1502 | // ============================================================================================ |
1503 | // the same for event mixing | |
1504 | AliTHn *fHistPMix = bfMix->GetHistNp(); | |
1505 | AliTHn *fHistNMix = bfMix->GetHistNn(); | |
1506 | AliTHn *fHistPNMix = bfMix->GetHistNpn(); | |
1507 | AliTHn *fHistNPMix = bfMix->GetHistNnp(); | |
1508 | AliTHn *fHistPPMix = bfMix->GetHistNpp(); | |
1509 | AliTHn *fHistNNMix = bfMix->GetHistNnn(); | |
1510 | ||
34616eda | 1511 | // pt trigger (fixed for all vertices/psi/centralities) |
f2e8af26 | 1512 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { |
622473b9 | 1513 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1514 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1515 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1516 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1517 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1518 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
f2e8af26 | 1519 | } |
1520 | ||
34616eda | 1521 | // pt associated (fixed for all vertices/psi/centralities) |
f2e8af26 | 1522 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) { |
622473b9 | 1523 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
1524 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1525 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
1526 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
f2e8af26 | 1527 | } |
f2e8af26 | 1528 | |
34616eda | 1529 | // ============================================================================================ |
1530 | // ranges (in bins) for vertexZ and centrality (psi) | |
1531 | Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin); | |
1532 | Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001); | |
a4e2c68f | 1533 | Int_t binVertexMin = 0; |
1534 | Int_t binVertexMax = 0; | |
1535 | if(fVertexBinning){ | |
1536 | binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin); | |
1537 | binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001); | |
1538 | } | |
34616eda | 1539 | Double_t binPsiLowEdge = 0.; |
1540 | Double_t binPsiUpEdge = 0.; | |
1541 | Double_t binVertexLowEdge = 0.; | |
1542 | Double_t binVertexUpEdge = 0.; | |
1543 | ||
1544 | TH2D* h1 = NULL; | |
1545 | TH2D* h2 = NULL; | |
1546 | TH2D* h3 = NULL; | |
1547 | TH2D* h4 = NULL; | |
1548 | ||
1549 | // loop over all bins | |
1550 | for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){ | |
1551 | for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){ | |
1552 | ||
1553 | cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl; | |
1554 | ||
1555 | // determine the bin edges for this bin | |
1556 | binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi); | |
1557 | binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi); | |
a4e2c68f | 1558 | if(fVertexBinning){ |
1559 | binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex); | |
1560 | binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex); | |
1561 | } | |
34616eda | 1562 | |
1563 | // Psi_2 | |
1564 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1565 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1566 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1567 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1568 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1569 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1570 | ||
1571 | // Vz | |
1572 | if(fVertexBinning){ | |
1573 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1574 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1575 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1576 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1577 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1578 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1579 | } | |
f2e8af26 | 1580 | |
34616eda | 1581 | // ============================================================================================ |
1582 | // the same for event mixing | |
1583 | ||
1584 | // Psi_2 | |
1585 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1586 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1587 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1588 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1589 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1590 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1591 | ||
1592 | // Vz | |
1593 | if(fVertexBinning){ | |
1594 | fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1595 | fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1596 | fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1597 | fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1598 | fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1599 | fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1600 | } | |
1601 | // ============================================================================================ | |
f2e8af26 | 1602 | |
34616eda | 1603 | //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0))); |
1604 | ||
1605 | // Project into the wanted space (1st: analysis step, 2nd: axis) | |
1606 | TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2); | |
1607 | TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2); | |
1608 | TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2); | |
1609 | TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2); | |
1610 | TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1); | |
1611 | TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1); | |
1612 | ||
1613 | // ============================================================================================ | |
1614 | // the same for event mixing | |
1615 | TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2); | |
1616 | TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2); | |
1617 | TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2); | |
1618 | TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2); | |
d631c24b | 1619 | // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1); |
1620 | // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1); | |
34616eda | 1621 | // ============================================================================================ |
1622 | ||
1623 | hTemp1->Sumw2(); | |
1624 | hTemp2->Sumw2(); | |
1625 | hTemp3->Sumw2(); | |
1626 | hTemp4->Sumw2(); | |
1627 | hTemp1Mix->Sumw2(); | |
1628 | hTemp2Mix->Sumw2(); | |
1629 | hTemp3Mix->Sumw2(); | |
1630 | hTemp4Mix->Sumw2(); | |
1631 | ||
1632 | hTemp1->Scale(1./hTemp5->GetEntries()); | |
1633 | hTemp3->Scale(1./hTemp5->GetEntries()); | |
1634 | hTemp2->Scale(1./hTemp6->GetEntries()); | |
1635 | hTemp4->Scale(1./hTemp6->GetEntries()); | |
271ceae3 | 1636 | |
1637 | // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method | |
1638 | Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX()); | |
1639 | mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1640 | hTemp1Mix->Scale(1./mixedNorm1); | |
1641 | Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX()); | |
1642 | mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1643 | hTemp2Mix->Scale(1./mixedNorm2); | |
1644 | Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX()); | |
1645 | mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1646 | hTemp3Mix->Scale(1./mixedNorm3); | |
1647 | Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX()); | |
1648 | mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1); | |
1649 | hTemp4Mix->Scale(1./mixedNorm4); | |
1650 | ||
34616eda | 1651 | hTemp1->Divide(hTemp1Mix); |
1652 | hTemp2->Divide(hTemp2Mix); | |
1653 | hTemp3->Divide(hTemp3Mix); | |
1654 | hTemp4->Divide(hTemp4Mix); | |
1655 | ||
1656 | // for the first: clone | |
1657 | if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){ | |
1658 | h1 = (TH2D*)hTemp1->Clone(); | |
1659 | h2 = (TH2D*)hTemp2->Clone(); | |
1660 | h3 = (TH2D*)hTemp3->Clone(); | |
1661 | h4 = (TH2D*)hTemp4->Clone(); | |
1662 | } | |
1663 | else{ // otherwise: add for averaging | |
1664 | h1->Add(hTemp1); | |
1665 | h2->Add(hTemp2); | |
1666 | h3->Add(hTemp3); | |
1667 | h4->Add(hTemp4); | |
1668 | } | |
f2e8af26 | 1669 | |
34616eda | 1670 | } |
1671 | } | |
f2e8af26 | 1672 | |
1673 | TH1D *gHistBalanceFunctionHistogram = 0x0; | |
34616eda | 1674 | if((h1)&&(h2)&&(h3)&&(h4)) { |
f2e8af26 | 1675 | |
34616eda | 1676 | // average over number of bins nbinsVertex * nbinsPsi |
1677 | h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1678 | h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1679 | h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
1680 | h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); | |
f2e8af26 | 1681 | |
1682 | // now only project on one axis | |
1683 | TH1D *h1DTemp1 = NULL; | |
1684 | TH1D *h1DTemp2 = NULL; | |
1685 | TH1D *h1DTemp3 = NULL; | |
1686 | TH1D *h1DTemp4 = NULL; | |
1687 | if(!bPhi){ | |
34616eda | 1688 | h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName())); |
1689 | h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName())); | |
1690 | h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName())); | |
1691 | h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName())); | |
f2e8af26 | 1692 | } |
1693 | else{ | |
34616eda | 1694 | h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName())); |
1695 | h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName())); | |
1696 | h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName())); | |
1697 | h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName())); | |
f2e8af26 | 1698 | } |
1699 | ||
1700 | gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName())); | |
1701 | gHistBalanceFunctionHistogram->Reset(); | |
1702 | if(!bPhi){ | |
1703 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta"); | |
1704 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)"); | |
1705 | } | |
1706 | else{ | |
1707 | gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)"); | |
1708 | gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)"); | |
1709 | } | |
1710 | ||
1711 | h1DTemp1->Add(h1DTemp3,-1.); | |
1712 | h1DTemp2->Add(h1DTemp4,-1.); | |
1713 | ||
1714 | gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.); | |
1715 | gHistBalanceFunctionHistogram->Scale(0.5); | |
a9f20288 | 1716 | |
1717 | //normalize to bin width | |
1718 | gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1))); | |
f2e8af26 | 1719 | } |
1720 | ||
1721 | return gHistBalanceFunctionHistogram; | |
1722 | } | |
1723 | ||
f0c5040c | 1724 | |
34616eda | 1725 | //____________________________________________________________________// |
1726 | TH2D *AliBalancePsi::GetCorrelationFunction(TString type, | |
1727 | Double_t psiMin, | |
1728 | Double_t psiMax, | |
1729 | Double_t vertexZMin, | |
1730 | Double_t vertexZMax, | |
1731 | Double_t ptTriggerMin, | |
1732 | Double_t ptTriggerMax, | |
1733 | Double_t ptAssociatedMin, | |
1734 | Double_t ptAssociatedMax, | |
1735 | AliBalancePsi *bMixed) { | |
1736 | ||
1737 | // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs, | |
1738 | // does the division by event mixing inside, | |
1739 | // and averaging over several vertexZ and centrality bins | |
1740 | ||
1741 | // security checks | |
1742 | if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){ | |
1743 | AliError("Only types allowed: PN,NP,PP,NN,ALL"); | |
1744 | return NULL; | |
1745 | } | |
1746 | if(!bMixed){ | |
1747 | AliError("No Event Mixing AliTHn"); | |
1748 | return NULL; | |
1749 | } | |
1750 | ||
1751 | TH2D *gHist = NULL; | |
1752 | TH2D *fSame = NULL; | |
1753 | TH2D *fMixed = NULL; | |
1754 | ||
1755 | // ranges (in bins) for vertexZ and centrality (psi) | |
1756 | Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin); | |
1757 | Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001); | |
a4e2c68f | 1758 | Int_t binVertexMin = 0; |
1759 | Int_t binVertexMax = 0; | |
1760 | if(fVertexBinning){ | |
1761 | binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin); | |
1762 | binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001); | |
1763 | } | |
34616eda | 1764 | Double_t binPsiLowEdge = 0.; |
a4e2c68f | 1765 | Double_t binPsiUpEdge = 1000.; |
34616eda | 1766 | Double_t binVertexLowEdge = 0.; |
a4e2c68f | 1767 | Double_t binVertexUpEdge = 1000.; |
34616eda | 1768 | |
1769 | // loop over all bins | |
1770 | for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){ | |
1771 | for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){ | |
1772 | ||
1773 | cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl; | |
1774 | ||
1775 | // determine the bin edges for this bin | |
1776 | binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi); | |
1777 | binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi); | |
a4e2c68f | 1778 | if(fVertexBinning){ |
1779 | binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex); | |
1780 | binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex); | |
1781 | } | |
1782 | ||
34616eda | 1783 | // get the 2D histograms for the correct type |
1784 | if(type=="PN"){ | |
1785 | fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1786 | fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1787 | } | |
1788 | else if(type=="NP"){ | |
1789 | fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1790 | fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1791 | } | |
1792 | else if(type=="PP"){ | |
1793 | fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
7071dbe3 | 1794 | fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); |
34616eda | 1795 | } |
1796 | else if(type=="NN"){ | |
1797 | fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1798 | fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1799 | } | |
1800 | else if(type=="ALL"){ | |
1801 | fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1802 | fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax); | |
1803 | } | |
1804 | ||
1fdb6a72 | 1805 | if(fSame && fMixed){ |
1806 | // then get the correlation function (divide fSame/fmixed) | |
1807 | fSame->Divide(fMixed); | |
1808 | ||
7071dbe3 | 1809 | // NEW averaging: |
1810 | // average over number of triggers in each sub-bin | |
a9737921 | 1811 | Double_t NTrigSubBin = 0; |
1812 | if(type=="PN" || type=="PP") | |
1813 | NTrigSubBin = (Double_t)(fHistP->Project(0,1)->GetEntries()); | |
1814 | else if(type=="NP" || type=="NN") | |
1815 | NTrigSubBin = (Double_t)(fHistN->Project(0,1)->GetEntries()); | |
7071dbe3 | 1816 | fSame->Scale(NTrigSubBin); |
1817 | ||
1fdb6a72 | 1818 | // for the first: clone |
84d14827 | 1819 | if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){ |
1fdb6a72 | 1820 | gHist = (TH2D*)fSame->Clone(); |
1821 | } | |
1822 | else{ // otherwise: add for averaging | |
1823 | gHist->Add(fSame); | |
1824 | } | |
34616eda | 1825 | } |
1826 | } | |
1827 | } | |
1828 | ||
1fdb6a72 | 1829 | if(gHist){ |
7071dbe3 | 1830 | |
1831 | // OLD averaging: | |
1fdb6a72 | 1832 | // average over number of bins nbinsVertex * nbinsPsi |
7071dbe3 | 1833 | // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1))); |
1834 | ||
1835 | // NEW averaging: | |
1836 | // average over number of triggers in each sub-bin | |
1837 | // first set to full range and then obtain number of all triggers | |
a9737921 | 1838 | Double_t NTrigAll = 0; |
1839 | if(type=="PN" || type=="PP"){ | |
1840 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1841 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1842 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1843 | NTrigAll = (Double_t)(fHistP->Project(0,1)->GetEntries()); | |
1844 | } | |
1845 | else if(type=="NP" || type=="NN"){ | |
1846 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
1847 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1848 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
1849 | NTrigAll = (Double_t)(fHistN->Project(0,1)->GetEntries()); | |
1850 | } | |
7071dbe3 | 1851 | gHist->Scale(1./NTrigAll); |
1852 | ||
1fdb6a72 | 1853 | } |
34616eda | 1854 | |
1855 | return gHist; | |
1856 | } | |
1857 | ||
1858 | ||
621329de | 1859 | //____________________________________________________________________// |
1860 | TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin, | |
6acdbcb2 | 1861 | Double_t psiMax, |
20006629 | 1862 | Double_t vertexZMin, |
1863 | Double_t vertexZMax, | |
6acdbcb2 | 1864 | Double_t ptTriggerMin, |
1865 | Double_t ptTriggerMax, | |
1866 | Double_t ptAssociatedMin, | |
47e040b5 | 1867 | Double_t ptAssociatedMax) { |
a38fd7f3 | 1868 | //Returns the 2D correlation function for (+-) pairs |
622473b9 | 1869 | |
1870 | // security checks | |
b0d53e32 | 1871 | if(psiMin > psiMax-0.00001){ |
1872 | AliError("psiMax <= psiMin"); | |
1873 | return NULL; | |
1874 | } | |
20006629 | 1875 | if(vertexZMin > vertexZMax-0.00001){ |
1876 | AliError("vertexZMax <= vertexZMin"); | |
1877 | return NULL; | |
1878 | } | |
b0d53e32 | 1879 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
1880 | AliError("ptTriggerMax <= ptTriggerMin"); | |
1881 | return NULL; | |
1882 | } | |
1883 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
1884 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
1885 | return NULL; | |
1886 | } | |
622473b9 | 1887 | |
621329de | 1888 | // Psi_2: axis 0 |
622473b9 | 1889 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
1890 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
621329de | 1891 | |
20006629 | 1892 | // Vz |
1893 | if(fVertexBinning){ | |
1894 | //Printf("Cutting on Vz..."); | |
1895 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1896 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1897 | } | |
1898 | ||
6acdbcb2 | 1899 | // pt trigger |
1900 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 1901 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1902 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
5ef8eaa7 | 1903 | } |
621329de | 1904 | |
6acdbcb2 | 1905 | // pt associated |
1906 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) | |
622473b9 | 1907 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
6acdbcb2 | 1908 | |
1909 | //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5); | |
1910 | //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5); | |
1911 | ||
1912 | //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1)); | |
1913 | //TCanvas *c1 = new TCanvas("c1",""); | |
1914 | //c1->cd(); | |
1915 | //if(!gHistTest){ | |
1916 | //AliError("Projection of fHistP = NULL"); | |
1917 | //return gHistTest; | |
1918 | //} | |
1919 | //else{ | |
1920 | //gHistTest->DrawCopy("colz"); | |
1921 | //} | |
1922 | ||
1923 | //0:step, 1: Delta eta, 2: Delta phi | |
621329de | 1924 | TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2)); |
89c00c43 | 1925 | if(!gHist){ |
1926 | AliError("Projection of fHistPN = NULL"); | |
1927 | return gHist; | |
1928 | } | |
1929 | ||
6acdbcb2 | 1930 | //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries()))); |
1931 | //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries()))); | |
1932 | //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries()))); | |
621329de | 1933 | |
6acdbcb2 | 1934 | //TCanvas *c2 = new TCanvas("c2",""); |
1935 | //c2->cd(); | |
1936 | //fHistPN->Project(0,1,2)->DrawCopy("colz"); | |
621329de | 1937 | |
47e040b5 | 1938 | if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0) |
621329de | 1939 | gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries())); |
faa03677 | 1940 | |
1941 | //normalize to bin width | |
1942 | gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1))); | |
a38fd7f3 | 1943 | |
1944 | return gHist; | |
1945 | } | |
1946 | ||
1947 | //____________________________________________________________________// | |
621329de | 1948 | TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin, |
6acdbcb2 | 1949 | Double_t psiMax, |
20006629 | 1950 | Double_t vertexZMin, |
1951 | Double_t vertexZMax, | |
6acdbcb2 | 1952 | Double_t ptTriggerMin, |
1953 | Double_t ptTriggerMax, | |
1954 | Double_t ptAssociatedMin, | |
47e040b5 | 1955 | Double_t ptAssociatedMax) { |
a38fd7f3 | 1956 | //Returns the 2D correlation function for (+-) pairs |
622473b9 | 1957 | |
1958 | // security checks | |
b0d53e32 | 1959 | if(psiMin > psiMax-0.00001){ |
1960 | AliError("psiMax <= psiMin"); | |
1961 | return NULL; | |
1962 | } | |
20006629 | 1963 | if(vertexZMin > vertexZMax-0.00001){ |
1964 | AliError("vertexZMax <= vertexZMin"); | |
1965 | return NULL; | |
1966 | } | |
b0d53e32 | 1967 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
1968 | AliError("ptTriggerMax <= ptTriggerMin"); | |
1969 | return NULL; | |
1970 | } | |
1971 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
1972 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
1973 | return NULL; | |
1974 | } | |
622473b9 | 1975 | |
621329de | 1976 | // Psi_2: axis 0 |
622473b9 | 1977 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
1978 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
20006629 | 1979 | |
1980 | // Vz | |
1981 | if(fVertexBinning){ | |
1982 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1983 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
1984 | } | |
a38fd7f3 | 1985 | |
6acdbcb2 | 1986 | // pt trigger |
1987 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 1988 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
1989 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
6acdbcb2 | 1990 | } |
1991 | ||
1992 | // pt associated | |
1993 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) | |
622473b9 | 1994 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
6acdbcb2 | 1995 | |
1996 | //0:step, 1: Delta eta, 2: Delta phi | |
621329de | 1997 | TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2)); |
89c00c43 | 1998 | if(!gHist){ |
1999 | AliError("Projection of fHistPN = NULL"); | |
2000 | return gHist; | |
2001 | } | |
2002 | ||
a38fd7f3 | 2003 | //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries())); |
2004 | //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries())); | |
47e040b5 | 2005 | if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0) |
621329de | 2006 | gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries())); |
faa03677 | 2007 | |
2008 | //normalize to bin width | |
2009 | gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1))); | |
a38fd7f3 | 2010 | |
2011 | return gHist; | |
2012 | } | |
2013 | ||
2014 | //____________________________________________________________________// | |
621329de | 2015 | TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin, |
6acdbcb2 | 2016 | Double_t psiMax, |
20006629 | 2017 | Double_t vertexZMin, |
2018 | Double_t vertexZMax, | |
6acdbcb2 | 2019 | Double_t ptTriggerMin, |
2020 | Double_t ptTriggerMax, | |
2021 | Double_t ptAssociatedMin, | |
47e040b5 | 2022 | Double_t ptAssociatedMax) { |
a38fd7f3 | 2023 | //Returns the 2D correlation function for (+-) pairs |
622473b9 | 2024 | |
2025 | // security checks | |
b0d53e32 | 2026 | if(psiMin > psiMax-0.00001){ |
2027 | AliError("psiMax <= psiMin"); | |
2028 | return NULL; | |
2029 | } | |
20006629 | 2030 | if(vertexZMin > vertexZMax-0.00001){ |
2031 | AliError("vertexZMax <= vertexZMin"); | |
2032 | return NULL; | |
2033 | } | |
b0d53e32 | 2034 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
2035 | AliError("ptTriggerMax <= ptTriggerMin"); | |
2036 | return NULL; | |
2037 | } | |
2038 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
2039 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
2040 | return NULL; | |
2041 | } | |
622473b9 | 2042 | |
621329de | 2043 | // Psi_2: axis 0 |
622473b9 | 2044 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
2045 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
6acdbcb2 | 2046 | |
20006629 | 2047 | // Vz |
2048 | if(fVertexBinning){ | |
2049 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2050 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2051 | } | |
2052 | ||
6acdbcb2 | 2053 | // pt trigger |
2054 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 2055 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
2056 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
6acdbcb2 | 2057 | } |
2058 | ||
2059 | // pt associated | |
2060 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) | |
622473b9 | 2061 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
a38fd7f3 | 2062 | |
6acdbcb2 | 2063 | //0:step, 1: Delta eta, 2: Delta phi |
621329de | 2064 | TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2)); |
89c00c43 | 2065 | if(!gHist){ |
2066 | AliError("Projection of fHistPN = NULL"); | |
2067 | return gHist; | |
2068 | } | |
2069 | ||
a38fd7f3 | 2070 | //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries())); |
2071 | //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries())); | |
47e040b5 | 2072 | if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0) |
621329de | 2073 | gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries())); |
faa03677 | 2074 | |
2075 | //normalize to bin width | |
2076 | gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1))); | |
a38fd7f3 | 2077 | |
2078 | return gHist; | |
2079 | } | |
2080 | ||
2081 | //____________________________________________________________________// | |
621329de | 2082 | TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin, |
6acdbcb2 | 2083 | Double_t psiMax, |
20006629 | 2084 | Double_t vertexZMin, |
2085 | Double_t vertexZMax, | |
6acdbcb2 | 2086 | Double_t ptTriggerMin, |
2087 | Double_t ptTriggerMax, | |
2088 | Double_t ptAssociatedMin, | |
47e040b5 | 2089 | Double_t ptAssociatedMax) { |
a38fd7f3 | 2090 | //Returns the 2D correlation function for (+-) pairs |
622473b9 | 2091 | |
2092 | // security checks | |
2093 | if(psiMin > psiMax-0.00001){ | |
2094 | AliError("psiMax <= psiMin"); | |
2095 | return NULL; | |
2096 | } | |
20006629 | 2097 | if(vertexZMin > vertexZMax-0.00001){ |
2098 | AliError("vertexZMax <= vertexZMin"); | |
2099 | return NULL; | |
2100 | } | |
622473b9 | 2101 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
2102 | AliError("ptTriggerMax <= ptTriggerMin"); | |
2103 | return NULL; | |
2104 | } | |
2105 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
2106 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
2107 | return NULL; | |
2108 | } | |
2109 | ||
621329de | 2110 | // Psi_2: axis 0 |
622473b9 | 2111 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
2112 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
6acdbcb2 | 2113 | |
20006629 | 2114 | // Vz |
2115 | if(fVertexBinning){ | |
2116 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2117 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2118 | } | |
2119 | ||
6acdbcb2 | 2120 | // pt trigger |
2121 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 2122 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
2123 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
6acdbcb2 | 2124 | } |
2125 | ||
2126 | // pt associated | |
2127 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) | |
622473b9 | 2128 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
a38fd7f3 | 2129 | |
6acdbcb2 | 2130 | //0:step, 1: Delta eta, 2: Delta phi |
621329de | 2131 | TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2)); |
89c00c43 | 2132 | if(!gHist){ |
2133 | AliError("Projection of fHistPN = NULL"); | |
2134 | return gHist; | |
2135 | } | |
2136 | ||
a38fd7f3 | 2137 | //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries())); |
2138 | //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries())); | |
47e040b5 | 2139 | if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0) |
621329de | 2140 | gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries())); |
faa03677 | 2141 | |
2142 | //normalize to bin width | |
2143 | gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1))); | |
a38fd7f3 | 2144 | |
2145 | return gHist; | |
2146 | } | |
621329de | 2147 | |
7b610f27 | 2148 | //____________________________________________________________________// |
2149 | TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin, | |
2150 | Double_t psiMax, | |
20006629 | 2151 | Double_t vertexZMin, |
2152 | Double_t vertexZMax, | |
7b610f27 | 2153 | Double_t ptTriggerMin, |
2154 | Double_t ptTriggerMax, | |
2155 | Double_t ptAssociatedMin, | |
47e040b5 | 2156 | Double_t ptAssociatedMax) { |
7b610f27 | 2157 | //Returns the 2D correlation function for the sum of all charge combination pairs |
622473b9 | 2158 | |
2159 | // security checks | |
b0d53e32 | 2160 | if(psiMin > psiMax-0.00001){ |
2161 | AliError("psiMax <= psiMin"); | |
2162 | return NULL; | |
2163 | } | |
20006629 | 2164 | if(vertexZMin > vertexZMax-0.00001){ |
2165 | AliError("vertexZMax <= vertexZMin"); | |
2166 | return NULL; | |
2167 | } | |
b0d53e32 | 2168 | if(ptTriggerMin > ptTriggerMax-0.00001){ |
2169 | AliError("ptTriggerMax <= ptTriggerMin"); | |
2170 | return NULL; | |
2171 | } | |
2172 | if(ptAssociatedMin > ptAssociatedMax-0.00001){ | |
2173 | AliError("ptAssociatedMax <= ptAssociatedMin"); | |
2174 | return NULL; | |
2175 | } | |
622473b9 | 2176 | |
7b610f27 | 2177 | // Psi_2: axis 0 |
622473b9 | 2178 | fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); |
2179 | fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
2180 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
2181 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
2182 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
2183 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001); | |
20006629 | 2184 | |
2185 | // Vz | |
2186 | if(fVertexBinning){ | |
2187 | fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2188 | fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2189 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2190 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2191 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2192 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001); | |
2193 | } | |
7b610f27 | 2194 | |
2195 | // pt trigger | |
2196 | if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) { | |
622473b9 | 2197 | fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); |
2198 | fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
2199 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
2200 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
2201 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
2202 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001); | |
7b610f27 | 2203 | } |
2204 | ||
2205 | // pt associated | |
2206 | if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) | |
622473b9 | 2207 | fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); |
2208 | fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
2209 | fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
2210 | fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001); | |
7b610f27 | 2211 | |
2212 | //0:step, 1: Delta eta, 2: Delta phi | |
2213 | TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2)); | |
2214 | if(!gHistNN){ | |
2215 | AliError("Projection of fHistNN = NULL"); | |
2216 | return gHistNN; | |
2217 | } | |
2218 | TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2)); | |
2219 | if(!gHistPP){ | |
2220 | AliError("Projection of fHistPP = NULL"); | |
2221 | return gHistPP; | |
2222 | } | |
2223 | TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2)); | |
2224 | if(!gHistNP){ | |
2225 | AliError("Projection of fHistNP = NULL"); | |
2226 | return gHistNP; | |
2227 | } | |
2228 | TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2)); | |
2229 | if(!gHistPN){ | |
2230 | AliError("Projection of fHistPN = NULL"); | |
2231 | return gHistPN; | |
2232 | } | |
2233 | ||
2234 | // sum all 2 particle histograms | |
2235 | gHistNN->Add(gHistPP); | |
2236 | gHistNN->Add(gHistNP); | |
2237 | gHistNN->Add(gHistPN); | |
2238 | ||
47e040b5 | 2239 | // divide by sum of + and - triggers |
2240 | if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0) | |
7b610f27 | 2241 | gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries())); |
faa03677 | 2242 | |
2243 | //normalize to bin width | |
2244 | gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1))); | |
7b610f27 | 2245 | |
2246 | return gHistNN; | |
2247 | } | |
2248 | ||
f2e8af26 | 2249 | //____________________________________________________________________// |
2250 | ||
1822002d | 2251 | Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM, |
f2e8af26 | 2252 | Double_t &mean, Double_t &meanError, |
2253 | Double_t &sigma, Double_t &sigmaError, | |
2254 | Double_t &skewness, Double_t &skewnessError, | |
2255 | Double_t &kurtosis, Double_t &kurtosisError) { | |
2256 | // | |
2257 | // helper method to calculate the moments and errors of a TH1D anlytically | |
4e023902 | 2258 | // fVariable = 1 for Delta eta (moments in full range) |
2259 | // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2) | |
f2e8af26 | 2260 | // |
2261 | ||
2262 | Bool_t success = kFALSE; | |
2263 | mean = -1.; | |
2264 | meanError = -1.; | |
2265 | sigma = -1.; | |
2266 | sigmaError = -1.; | |
2267 | skewness = -1.; | |
2268 | skewnessError = -1.; | |
2269 | kurtosis = -1.; | |
2270 | kurtosisError = -1.; | |
2271 | ||
2272 | if(gHist){ | |
2273 | ||
2274 | // ---------------------------------------------------------------------- | |
2275 | // basic parameters of histogram | |
2276 | ||
2277 | Int_t fNumberOfBins = gHist->GetNbinsX(); | |
2278 | // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning | |
2279 | ||
1a879e25 | 2280 | |
2281 | // ---------------------------------------------------------------------- | |
2282 | // ZYAM (for partially negative distributions) | |
2283 | // --> we subtract always the minimum value | |
7071dbe3 | 2284 | Double_t zeroYield = 0.; |
2285 | Double_t zeroYieldCur = -FLT_MAX; | |
2286 | if(kUseZYAM){ | |
4184be06 | 2287 | for(Int_t iMin = 0; iMin<2; iMin++){ |
7071dbe3 | 2288 | zeroYieldCur = gHist->GetMinimum(zeroYieldCur); |
2289 | zeroYield += zeroYieldCur; | |
2290 | } | |
4184be06 | 2291 | zeroYield /= 2.; |
7071dbe3 | 2292 | //zeroYield = gHist->GetMinimum(); |
2293 | } | |
f2e8af26 | 2294 | // ---------------------------------------------------------------------- |
2295 | // first calculate the mean | |
2296 | ||
2297 | Double_t fWeightedAverage = 0.; | |
f0754617 | 2298 | Double_t fNormalization = 0.; |
f2e8af26 | 2299 | |
2300 | for(Int_t i = 1; i <= fNumberOfBins; i++) { | |
2301 | ||
4e023902 | 2302 | // for Delta phi: moments only on near side -pi/2 < dphi < pi/2 |
2303 | if(fVariable == 2 && | |
2304 | (gHist->GetBinCenter(i) < - TMath::Pi()/2 || | |
2305 | gHist->GetBinCenter(i) > TMath::Pi()/2)){ | |
2306 | continue; | |
2307 | } | |
2308 | ||
1a879e25 | 2309 | fWeightedAverage += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i); |
2310 | fNormalization += (gHist->GetBinContent(i)-zeroYield); | |
f2e8af26 | 2311 | } |
2312 | ||
2313 | mean = fWeightedAverage / fNormalization; | |
2314 | ||
f2e8af26 | 2315 | // ---------------------------------------------------------------------- |
2316 | // then calculate the higher moments | |
2317 | ||
a123fe62 | 2318 | Double_t fMu = 0.; |
2319 | Double_t fMu2 = 0.; | |
2320 | Double_t fMu3 = 0.; | |
2321 | Double_t fMu4 = 0.; | |
2322 | Double_t fMu5 = 0.; | |
2323 | Double_t fMu6 = 0.; | |
2324 | Double_t fMu7 = 0.; | |
2325 | Double_t fMu8 = 0.; | |
f2e8af26 | 2326 | |
2327 | for(Int_t i = 1; i <= fNumberOfBins; i++) { | |
4e023902 | 2328 | |
2329 | // for Delta phi: moments only on near side -pi/2 < dphi < pi/2 | |
2330 | if(fVariable == 2 && | |
2331 | (gHist->GetBinCenter(i) < - TMath::Pi()/2 || | |
2332 | gHist->GetBinCenter(i) > TMath::Pi()/2)){ | |
2333 | continue; | |
2334 | } | |
2335 | ||
1a879e25 | 2336 | fMu += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean); |
2337 | fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2); | |
2338 | fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3); | |
2339 | fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4); | |
2340 | fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5); | |
2341 | fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6); | |
2342 | fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7); | |
2343 | fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8); | |
f2e8af26 | 2344 | } |
e47eeabd | 2345 | |
2346 | // normalize to bin entries! | |
2347 | fMu /= fNormalization; | |
2348 | fMu2 /= fNormalization; | |
2349 | fMu3 /= fNormalization; | |
2350 | fMu4 /= fNormalization; | |
2351 | fMu5 /= fNormalization; | |
2352 | fMu6 /= fNormalization; | |
2353 | fMu7 /= fNormalization; | |
2354 | fMu8 /= fNormalization; | |
2355 | ||
2356 | sigma = TMath::Sqrt(fMu2); | |
2357 | skewness = fMu3 / TMath::Power(sigma,3); | |
2358 | kurtosis = fMu4 / TMath::Power(sigma,4) - 3; | |
f2e8af26 | 2359 | |
efef044a | 2360 | // normalize with sigma^r |
2361 | fMu /= TMath::Power(sigma,1); | |
2362 | fMu2 /= TMath::Power(sigma,2); | |
2363 | fMu3 /= TMath::Power(sigma,3); | |
2364 | fMu4 /= TMath::Power(sigma,4); | |
2365 | fMu5 /= TMath::Power(sigma,5); | |
2366 | fMu6 /= TMath::Power(sigma,6); | |
2367 | fMu7 /= TMath::Power(sigma,7); | |
2368 | fMu8 /= TMath::Power(sigma,8); | |
1a879e25 | 2369 | |
f0754617 | 2370 | // ---------------------------------------------------------------------- |
a123fe62 | 2371 | // then calculate the higher moment errors |
4e023902 | 2372 | // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl; |
2373 | // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl; | |
2374 | // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl; | |
a123fe62 | 2375 | |
2376 | Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError | |
2377 | ||
2378 | // // assuming normal distribution (as done in ROOT) | |
2379 | // meanError = sigma / TMath::Sqrt(normError); | |
2380 | // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError)); | |
2381 | // skewnessError = TMath::Sqrt(6./(normError)); | |
2382 | // kurtosisError = TMath::Sqrt(24./(normError)); | |
2383 | ||
2384 | // use delta theorem paper (Luo - arXiv:1109.0593v1) | |
a9f20288 | 2385 | Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError)); |
2386 | Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError); | |
2387 | Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError); | |
1af1cefb | 2388 | //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError); |
2389 | //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError); | |
2390 | //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError; | |
a123fe62 | 2391 | |
4e023902 | 2392 | // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl; |
2393 | // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl; | |
a123fe62 | 2394 | |
a9f20288 | 2395 | if (TMath::Sqrt(normError) != 0){ |
2396 | meanError = sigma / TMath::Sqrt(normError); | |
1a879e25 | 2397 | sigmaError = TMath::Sqrt(Lambda11); |
2398 | skewnessError = TMath::Sqrt(Lambda22); | |
2399 | kurtosisError = TMath::Sqrt(Lambda33); | |
2400 | ||
2401 | success = kTRUE; | |
2402 | ||
a9f20288 | 2403 | } |
1a879e25 | 2404 | else success = kFALSE; |
2405 | ||
f2e8af26 | 2406 | } |
f2e8af26 | 2407 | return success; |
2408 | } | |
2409 | ||
2410 | ||
73609a48 | 2411 | //____________________________________________________________________// |
2412 | Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) { | |
2413 | // | |
2414 | // calculates dphistar | |
2415 | // | |
2416 | Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2); | |
2417 | ||
2418 | static const Double_t kPi = TMath::Pi(); | |
2419 | ||
2420 | // circularity | |
2421 | // if (dphistar > 2 * kPi) | |
2422 | // dphistar -= 2 * kPi; | |
2423 | // if (dphistar < -2 * kPi) | |
2424 | // dphistar += 2 * kPi; | |
2425 | ||
2426 | if (dphistar > kPi) | |
2427 | dphistar = kPi * 2 - dphistar; | |
2428 | if (dphistar < -kPi) | |
2429 | dphistar = -kPi * 2 - dphistar; | |
2430 | if (dphistar > kPi) // might look funny but is needed | |
2431 | dphistar = kPi * 2 - dphistar; | |
2432 | ||
2433 | return dphistar; | |
2434 | } | |
2435 | ||
2922bbe3 | 2436 | |
2437 |