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