]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coverity fixes - Part I (compilation warnings: set but not used)
[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
f0c5040c 797//____________________________________________________________________//
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
a38fd7f3 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
f0c5040c 1340//____________________________________________________________________//
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
1700 TH2D* h1 = NULL;
1701 TH2D* h2 = NULL;
1702 TH2D* h3 = NULL;
1703 TH2D* h4 = NULL;
1704
1705 // loop over all bins
1706 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1707 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1708
1709 cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1710
34616eda 1711 // Psi_2
1712 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1713 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1714 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1715 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1716 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1717 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1718
1719 // Vz
1720 if(fVertexBinning){
1721 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1722 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1723 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1724 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1725 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1726 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1727 }
f2e8af26 1728
34616eda 1729 // ============================================================================================
1730 // the same for event mixing
1731
1732 // Psi_2
1733 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1734 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1735 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1736 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1737 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1738 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1739
1740 // Vz
1741 if(fVertexBinning){
1742 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1743 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1744 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1745 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1746 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1747 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1748 }
1749 // ============================================================================================
f2e8af26 1750
34616eda 1751 //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)));
1752
1753 // Project into the wanted space (1st: analysis step, 2nd: axis)
1754 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1755 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1756 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1757 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1758 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1759 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1760
1761 // ============================================================================================
1762 // the same for event mixing
1763 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1764 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1765 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1766 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
d631c24b 1767 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1768 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
34616eda 1769 // ============================================================================================
1770
1771 hTemp1->Sumw2();
1772 hTemp2->Sumw2();
1773 hTemp3->Sumw2();
1774 hTemp4->Sumw2();
1775 hTemp1Mix->Sumw2();
1776 hTemp2Mix->Sumw2();
1777 hTemp3Mix->Sumw2();
1778 hTemp4Mix->Sumw2();
1779
5bb3cd60
MW		 1780 hTemp1->Scale(1./hTemp5->Integral());
1781 hTemp3->Scale(1./hTemp5->Integral());
1782 hTemp2->Scale(1./hTemp6->Integral());
1783 hTemp4->Scale(1./hTemp6->Integral());
271ceae3 1784
1785 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1786 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1787 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1788 hTemp1Mix->Scale(1./mixedNorm1);
1789 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1790 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1791 hTemp2Mix->Scale(1./mixedNorm2);
1792 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1793 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1794 hTemp3Mix->Scale(1./mixedNorm3);
1795 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1796 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1797 hTemp4Mix->Scale(1./mixedNorm4);
1798
34616eda 1799 hTemp1->Divide(hTemp1Mix);
1800 hTemp2->Divide(hTemp2Mix);
1801 hTemp3->Divide(hTemp3Mix);
1802 hTemp4->Divide(hTemp4Mix);
1803
1804 // for the first: clone
1805 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1806 h1 = (TH2D*)hTemp1->Clone();
1807 h2 = (TH2D*)hTemp2->Clone();
1808 h3 = (TH2D*)hTemp3->Clone();
1809 h4 = (TH2D*)hTemp4->Clone();
1810 }
1811 else{ // otherwise: add for averaging
1812 h1->Add(hTemp1);
1813 h2->Add(hTemp2);
1814 h3->Add(hTemp3);
1815 h4->Add(hTemp4);
1816 }
f2e8af26 1817
34616eda 1818 }
1819 }
f2e8af26 1820
1821 TH1D *gHistBalanceFunctionHistogram = 0x0;
34616eda 1822 if((h1)&&(h2)&&(h3)&&(h4)) {
f2e8af26 1823
34616eda 1824 // average over number of bins nbinsVertex * nbinsPsi
1825 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1826 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1827 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1828 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
f2e8af26 1829
1830 // now only project on one axis
1831 TH1D *h1DTemp1 = NULL;
1832 TH1D *h1DTemp2 = NULL;
1833 TH1D *h1DTemp3 = NULL;
1834 TH1D *h1DTemp4 = NULL;
1835 if(!bPhi){
34616eda 1836 h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName()));
1837 h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName()));
1838 h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName()));
1839 h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName()));
f2e8af26 1840 }
1841 else{
34616eda 1842 h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName()));
1843 h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName()));
1844 h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName()));
1845 h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName()));
f2e8af26 1846 }
1847
1848 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1849 gHistBalanceFunctionHistogram->Reset();
1850 if(!bPhi){
1851 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1852 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1853 }
1854 else{
1855 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1856 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1857 }
1858
1859 h1DTemp1->Add(h1DTemp3,-1.);
1860 h1DTemp2->Add(h1DTemp4,-1.);
1861
1862 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1863 gHistBalanceFunctionHistogram->Scale(0.5);
a9f20288 1864
1865 //normalize to bin width
1866 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
f2e8af26 1867 }
1868
1869 return gHistBalanceFunctionHistogram;
1870}
1871
f0c5040c 1872
34616eda 1873//____________________________________________________________________//
1874TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1875 Double_t psiMin,
1876 Double_t psiMax,
1877 Double_t vertexZMin,
1878 Double_t vertexZMax,
1879 Double_t ptTriggerMin,
1880 Double_t ptTriggerMax,
1881 Double_t ptAssociatedMin,
1882 Double_t ptAssociatedMax,
bd7bf11d 1883 AliBalancePsi *bMixed,
94415ad1 1884 Bool_t normToTrig,
c965fa95 1885 Double_t normalizationRangePhi) {
34616eda 1886
1887 // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1888 // does the division by event mixing inside,
1889 // and averaging over several vertexZ and centrality bins
1890
1891 // security checks
1892 if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1893 AliError("Only types allowed: PN,NP,PP,NN,ALL");
1894 return NULL;
1895 }
1896 if(!bMixed){
1897 AliError("No Event Mixing AliTHn");
1898 return NULL;
1899 }
1900
1901 TH2D *gHist = NULL;
1902 TH2D *fSame = NULL;
1903 TH2D *fMixed = NULL;
1904
1905 // ranges (in bins) for vertexZ and centrality (psi)
94415ad1 1906 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
34616eda 1907 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
a4e2c68f 1908 Int_t binVertexMin = 0;
1909 Int_t binVertexMax = 0;
1910 if(fVertexBinning){
94415ad1 1911 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
a4e2c68f 1912 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1913 }
34616eda 1914 Double_t binPsiLowEdge = 0.;
a4e2c68f 1915 Double_t binPsiUpEdge = 1000.;
34616eda 1916 Double_t binVertexLowEdge = 0.;
a4e2c68f 1917 Double_t binVertexUpEdge = 1000.;
34616eda 1918
1919 // loop over all bins
1920 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1921 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1922
3aa5487c 1923 //cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
34616eda 1924
1925 // determine the bin edges for this bin
94415ad1 1926 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
1927 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
a4e2c68f 1928 if(fVertexBinning){
94415ad1 1929 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
1930 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
a4e2c68f 1931 }
1932
34616eda 1933 // get the 2D histograms for the correct type
1934 if(type=="PN"){
1935 fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1936 fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1937 }
1938 else if(type=="NP"){
1939 fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1940 fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1941 }
1942 else if(type=="PP"){
1943 fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
7071dbe3 1944 fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
34616eda 1945 }
1946 else if(type=="NN"){
1947 fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1948 fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1949 }
1950 else if(type=="ALL"){
1951 fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1952 fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1953 }
94415ad1 1954
f7e8f569 1955 if(fMixed && normToTrig && fMixed->Integral()>0){
bd7bf11d
MW		 1956
1957 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
3aa5487c 1958 // do it only on away-side (due to two-track cuts): pi +- pi/6.
51bfb3bc 1959 Int_t binXmin = fMixed->GetXaxis()->FindBin(0-10e-5);
1960 Int_t binXmax = fMixed->GetXaxis()->FindBin(0+10e-5);
1961 Double_t binsX = (Double_t)(binXmax - binXmin + 1);
3aa5487c 1962 Int_t binYmin = fMixed->GetYaxis()->FindBin(TMath::Pi() - normalizationRangePhi);
1963 Int_t binYmax = fMixed->GetYaxis()->FindBin(TMath::Pi() + normalizationRangePhi - 0.00001);
51bfb3bc 1964 Double_t binsY = (Double_t)(binYmax - binYmin + 1);
1965
1966 Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
1967 mixedNorm /= binsX * binsY;
bd7bf11d
MW		 1968
1969 // finite bin correction
1970 Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
1971 Double_t maxEta = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
1972
1973 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
5bb3cd60 1974 //Printf("Finite bin correction: %f", finiteBinCorrection);
bd7bf11d
MW		 1975 mixedNorm /= finiteBinCorrection;
1976
1977 fMixed->Scale(1./mixedNorm);
1978 }
1979
1fdb6a72 1980 if(fSame && fMixed){
1981 // then get the correlation function (divide fSame/fmixed)
1982 fSame->Divide(fMixed);
1983
c965fa95 1984 // averaging with number of triggers:
1985 // average over number of triggers in each sub-bin
1986 Double_t NTrigSubBin = 0;
1987 if(type=="PN" || type=="PP")
1988 NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
1989 else if(type=="NP" || type=="NN")
1990 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
e6023587 1991 else if(type=="ALL")
1992 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
c965fa95 1993 fSame->Scale(NTrigSubBin);
7071dbe3 1994
1fdb6a72 1995 // for the first: clone
84d14827 1996 if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
1fdb6a72 1997 gHist = (TH2D*)fSame->Clone();
1998 }
1999 else{ // otherwise: add for averaging
2000 gHist->Add(fSame);
2001 }
34616eda 2002 }
2003 }
2004 }
2005
1fdb6a72 2006 if(gHist){
7071dbe3 2007
c965fa95 2008 // averaging with number of triggers:
2009 // first set to full range and then obtain number of all triggers
2010 Double_t NTrigAll = 0;
2011 if(type=="PN" || type=="PP"){
2012 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2013 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2014 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2015 NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
2016 }
2017 else if(type=="NP" || type=="NN"){
2018 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2019 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2020 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2021 NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
2022 }
e6023587 2023 else if(type=="ALL"){
2024 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2025 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2026 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2027 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2028 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2029 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2030 NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
2031 }
c965fa95 2032 gHist->Scale(1./NTrigAll);
2033
1fdb6a72 2034 }
34616eda 2035
2036 return gHist;
2037}
2038
2039
621329de 2040//____________________________________________________________________//
2041TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 2042 Double_t psiMax,
20006629 2043 Double_t vertexZMin,
2044 Double_t vertexZMax,
6acdbcb2 2045 Double_t ptTriggerMin,
2046 Double_t ptTriggerMax,
2047 Double_t ptAssociatedMin,
47e040b5 2048 Double_t ptAssociatedMax) {
a38fd7f3 2049 //Returns the 2D correlation function for (+-) pairs
622473b9 2050
2051 // security checks
b0d53e32 2052 if(psiMin > psiMax-0.00001){
2053 AliError("psiMax <= psiMin");
2054 return NULL;
2055 }
20006629 2056 if(vertexZMin > vertexZMax-0.00001){
2057 AliError("vertexZMax <= vertexZMin");
2058 return NULL;
2059 }
b0d53e32 2060 if(ptTriggerMin > ptTriggerMax-0.00001){
2061 AliError("ptTriggerMax <= ptTriggerMin");
2062 return NULL;
2063 }
2064 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2065 AliError("ptAssociatedMax <= ptAssociatedMin");
2066 return NULL;
2067 }
622473b9 2068
621329de 2069 // Psi_2: axis 0
622473b9 2070 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2071 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 2072
20006629 2073 // Vz
2074 if(fVertexBinning){
2075 //Printf("Cutting on Vz...");
2076 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2077 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2078 }
2079
6acdbcb2 2080 // pt trigger
2081 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2082 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2083 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 2084 }
621329de 2085
6acdbcb2 2086 // pt associated
2087 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2088 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 2089
2090 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2091 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2092
2093 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
2094 //TCanvas *c1 = new TCanvas("c1","");
2095 //c1->cd();
2096 //if(!gHistTest){
2097 //AliError("Projection of fHistP = NULL");
2098 //return gHistTest;
2099 //}
2100 //else{
2101 //gHistTest->DrawCopy("colz");
2102 //}
2103
2104 //0:step, 1: Delta eta, 2: Delta phi
621329de 2105 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 2106 if(!gHist){
2107 AliError("Projection of fHistPN = NULL");
2108 return gHist;
2109 }
2110
6acdbcb2 2111 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
2112 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
2113 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
5bb3cd60
MW		 2114 //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
2115 //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
621329de 2116
6acdbcb2 2117 //TCanvas *c2 = new TCanvas("c2","");
2118 //c2->cd();
2119 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 2120
5bb3cd60
MW		 2121 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2122 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
faa03677 2123
2124 //normalize to bin width
2125 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2126
2127 return gHist;
2128}
2129
2130//____________________________________________________________________//
621329de 2131TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 2132 Double_t psiMax,
20006629 2133 Double_t vertexZMin,
2134 Double_t vertexZMax,
6acdbcb2 2135 Double_t ptTriggerMin,
2136 Double_t ptTriggerMax,
2137 Double_t ptAssociatedMin,
47e040b5 2138 Double_t ptAssociatedMax) {
a38fd7f3 2139 //Returns the 2D correlation function for (+-) pairs
622473b9 2140
2141 // security checks
b0d53e32 2142 if(psiMin > psiMax-0.00001){
2143 AliError("psiMax <= psiMin");
2144 return NULL;
2145 }
20006629 2146 if(vertexZMin > vertexZMax-0.00001){
2147 AliError("vertexZMax <= vertexZMin");
2148 return NULL;
2149 }
b0d53e32 2150 if(ptTriggerMin > ptTriggerMax-0.00001){
2151 AliError("ptTriggerMax <= ptTriggerMin");
2152 return NULL;
2153 }
2154 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2155 AliError("ptAssociatedMax <= ptAssociatedMin");
2156 return NULL;
2157 }
622473b9 2158
621329de 2159 // Psi_2: axis 0
622473b9 2160 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2161 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 2162
2163 // Vz
2164 if(fVertexBinning){
2165 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2166 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2167 }
a38fd7f3 2168
6acdbcb2 2169 // pt trigger
2170 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2171 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2172 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2173 }
2174
2175 // pt associated
2176 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2177 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 2178
2179 //0:step, 1: Delta eta, 2: Delta phi
621329de 2180 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 2181 if(!gHist){
2182 AliError("Projection of fHistPN = NULL");
2183 return gHist;
2184 }
2185
a38fd7f3 2186 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2187 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
5bb3cd60
MW		 2188 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2189 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
faa03677 2190
2191 //normalize to bin width
2192 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2193
2194 return gHist;
2195}
2196
2197//____________________________________________________________________//
621329de 2198TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 2199 Double_t psiMax,
20006629 2200 Double_t vertexZMin,
2201 Double_t vertexZMax,
6acdbcb2 2202 Double_t ptTriggerMin,
2203 Double_t ptTriggerMax,
2204 Double_t ptAssociatedMin,
47e040b5 2205 Double_t ptAssociatedMax) {
a38fd7f3 2206 //Returns the 2D correlation function for (+-) pairs
622473b9 2207
2208 // security checks
b0d53e32 2209 if(psiMin > psiMax-0.00001){
2210 AliError("psiMax <= psiMin");
2211 return NULL;
2212 }
20006629 2213 if(vertexZMin > vertexZMax-0.00001){
2214 AliError("vertexZMax <= vertexZMin");
2215 return NULL;
2216 }
b0d53e32 2217 if(ptTriggerMin > ptTriggerMax-0.00001){
2218 AliError("ptTriggerMax <= ptTriggerMin");
2219 return NULL;
2220 }
2221 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2222 AliError("ptAssociatedMax <= ptAssociatedMin");
2223 return NULL;
2224 }
622473b9 2225
621329de 2226 // Psi_2: axis 0
622473b9 2227 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2228 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 2229
20006629 2230 // Vz
2231 if(fVertexBinning){
2232 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2233 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2234 }
2235
6acdbcb2 2236 // pt trigger
2237 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2238 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2239 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2240 }
2241
2242 // pt associated
2243 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2244 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 2245
6acdbcb2 2246 //0:step, 1: Delta eta, 2: Delta phi
621329de 2247 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 2248 if(!gHist){
2249 AliError("Projection of fHistPN = NULL");
2250 return gHist;
2251 }
2252
a38fd7f3 2253 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2254 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
5bb3cd60
MW		 2255 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2256 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
faa03677 2257
2258 //normalize to bin width
2259 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2260
2261 return gHist;
2262}
2263
2264//____________________________________________________________________//
621329de 2265TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 2266 Double_t psiMax,
20006629 2267 Double_t vertexZMin,
2268 Double_t vertexZMax,
6acdbcb2 2269 Double_t ptTriggerMin,
2270 Double_t ptTriggerMax,
2271 Double_t ptAssociatedMin,
47e040b5 2272 Double_t ptAssociatedMax) {
a38fd7f3 2273 //Returns the 2D correlation function for (+-) pairs
622473b9 2274
2275 // security checks
2276 if(psiMin > psiMax-0.00001){
2277 AliError("psiMax <= psiMin");
2278 return NULL;
2279 }
20006629 2280 if(vertexZMin > vertexZMax-0.00001){
2281 AliError("vertexZMax <= vertexZMin");
2282 return NULL;
2283 }
622473b9 2284 if(ptTriggerMin > ptTriggerMax-0.00001){
2285 AliError("ptTriggerMax <= ptTriggerMin");
2286 return NULL;
2287 }
2288 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2289 AliError("ptAssociatedMax <= ptAssociatedMin");
2290 return NULL;
2291 }
2292
621329de 2293 // Psi_2: axis 0
622473b9 2294 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2295 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 2296
20006629 2297 // Vz
2298 if(fVertexBinning){
2299 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2300 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2301 }
2302
6acdbcb2 2303 // pt trigger
2304 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2305 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2306 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2307 }
2308
2309 // pt associated
2310 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2311 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 2312
6acdbcb2 2313 //0:step, 1: Delta eta, 2: Delta phi
621329de 2314 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 2315 if(!gHist){
2316 AliError("Projection of fHistPN = NULL");
2317 return gHist;
2318 }
2319
a38fd7f3 2320 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2321 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
5bb3cd60
MW		 2322 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2323 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
faa03677 2324
2325 //normalize to bin width
2326 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2327
2328 return gHist;
2329}
621329de 2330
7b610f27 2331//____________________________________________________________________//
2332TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2333 Double_t psiMax,
20006629 2334 Double_t vertexZMin,
2335 Double_t vertexZMax,
7b610f27 2336 Double_t ptTriggerMin,
2337 Double_t ptTriggerMax,
2338 Double_t ptAssociatedMin,
47e040b5 2339 Double_t ptAssociatedMax) {
7b610f27 2340 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 2341
2342 // security checks
b0d53e32 2343 if(psiMin > psiMax-0.00001){
2344 AliError("psiMax <= psiMin");
2345 return NULL;
2346 }
20006629 2347 if(vertexZMin > vertexZMax-0.00001){
2348 AliError("vertexZMax <= vertexZMin");
2349 return NULL;
2350 }
b0d53e32 2351 if(ptTriggerMin > ptTriggerMax-0.00001){
2352 AliError("ptTriggerMax <= ptTriggerMin");
2353 return NULL;
2354 }
2355 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2356 AliError("ptAssociatedMax <= ptAssociatedMin");
2357 return NULL;
2358 }
622473b9 2359
7b610f27 2360 // Psi_2: axis 0
622473b9 2361 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2362 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2363 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2364 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2365 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2366 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 2367
2368 // Vz
2369 if(fVertexBinning){
2370 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2371 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2372 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2373 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2374 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2375 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2376 }
7b610f27 2377
2378 // pt trigger
2379 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2380 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2381 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2382 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2383 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2384 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2385 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 2386 }
2387
2388 // pt associated
2389 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2390 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2391 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2392 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2393 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 2394
2395 //0:step, 1: Delta eta, 2: Delta phi
2396 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2397 if(!gHistNN){
2398 AliError("Projection of fHistNN = NULL");
2399 return gHistNN;
2400 }
2401 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2402 if(!gHistPP){
2403 AliError("Projection of fHistPP = NULL");
2404 return gHistPP;
2405 }
2406 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2407 if(!gHistNP){
2408 AliError("Projection of fHistNP = NULL");
2409 return gHistNP;
2410 }
2411 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2412 if(!gHistPN){
2413 AliError("Projection of fHistPN = NULL");
2414 return gHistPN;
2415 }
2416
2417 // sum all 2 particle histograms
2418 gHistNN->Add(gHistPP);
2419 gHistNN->Add(gHistNP);
2420 gHistNN->Add(gHistPN);
2421
47e040b5 2422 // divide by sum of + and - triggers
5bb3cd60 2423 if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
e6023587 2424 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral()));
faa03677 2425
2426 //normalize to bin width
2427 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 2428
2429 return gHistNN;
2430}
2431
f2e8af26 2432//____________________________________________________________________//
2433
1822002d 2434Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
f2e8af26 2435 Double_t &mean, Double_t &meanError,
2436 Double_t &sigma, Double_t &sigmaError,
2437 Double_t &skewness, Double_t &skewnessError,
2438 Double_t &kurtosis, Double_t &kurtosisError) {
2439 //
2440 // helper method to calculate the moments and errors of a TH1D anlytically
4e023902 2441 // fVariable = 1 for Delta eta (moments in full range)
2442 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
f2e8af26 2443 //
2444
2445 Bool_t success = kFALSE;
2446 mean = -1.;
2447 meanError = -1.;
2448 sigma = -1.;
2449 sigmaError = -1.;
2450 skewness = -1.;
2451 skewnessError = -1.;
2452 kurtosis = -1.;
2453 kurtosisError = -1.;
2454
2455 if(gHist){
2456
2457 // ----------------------------------------------------------------------
2458 // basic parameters of histogram
2459
2460 Int_t fNumberOfBins = gHist->GetNbinsX();
2461 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
2462
1a879e25 2463
2464 // ----------------------------------------------------------------------
2465 // ZYAM (for partially negative distributions)
2466 // --> we subtract always the minimum value
7071dbe3 2467 Double_t zeroYield = 0.;
2468 Double_t zeroYieldCur = -FLT_MAX;
2469 if(kUseZYAM){
4184be06 2470 for(Int_t iMin = 0; iMin<2; iMin++){
7071dbe3 2471 zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2472 zeroYield += zeroYieldCur;
2473 }
4184be06 2474 zeroYield /= 2.;
7071dbe3 2475 //zeroYield = gHist->GetMinimum();
2476 }
f2e8af26 2477 // ----------------------------------------------------------------------
2478 // first calculate the mean
2479
2480 Double_t fWeightedAverage = 0.;
f0754617 2481 Double_t fNormalization = 0.;
f2e8af26 2482
2483 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2484
4e023902 2485 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2486 if(fVariable == 2 &&
2487 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2488 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2489 continue;
2490 }
2491
1a879e25 2492 fWeightedAverage += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2493 fNormalization += (gHist->GetBinContent(i)-zeroYield);
f2e8af26 2494 }
2495
2496 mean = fWeightedAverage / fNormalization;
2497
f2e8af26 2498 // ----------------------------------------------------------------------
2499 // then calculate the higher moments
2500
a123fe62 2501 Double_t fMu = 0.;
2502 Double_t fMu2 = 0.;
2503 Double_t fMu3 = 0.;
2504 Double_t fMu4 = 0.;
2505 Double_t fMu5 = 0.;
2506 Double_t fMu6 = 0.;
2507 Double_t fMu7 = 0.;
2508 Double_t fMu8 = 0.;
f2e8af26 2509
2510 for(Int_t i = 1; i <= fNumberOfBins; i++) {
4e023902 2511
2512 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2513 if(fVariable == 2 &&
2514 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2515 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2516 continue;
2517 }
2518
1a879e25 2519 fMu += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2520 fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2521 fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2522 fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2523 fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2524 fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2525 fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2526 fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
f2e8af26 2527 }
e47eeabd 2528
2529 // normalize to bin entries!
2530 fMu /= fNormalization;
2531 fMu2 /= fNormalization;
2532 fMu3 /= fNormalization;
2533 fMu4 /= fNormalization;
2534 fMu5 /= fNormalization;
2535 fMu6 /= fNormalization;
2536 fMu7 /= fNormalization;
2537 fMu8 /= fNormalization;
2538
2539 sigma = TMath::Sqrt(fMu2);
2540 skewness = fMu3 / TMath::Power(sigma,3);
2541 kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
f2e8af26 2542
efef044a 2543 // normalize with sigma^r
2544 fMu /= TMath::Power(sigma,1);
2545 fMu2 /= TMath::Power(sigma,2);
2546 fMu3 /= TMath::Power(sigma,3);
2547 fMu4 /= TMath::Power(sigma,4);
2548 fMu5 /= TMath::Power(sigma,5);
2549 fMu6 /= TMath::Power(sigma,6);
2550 fMu7 /= TMath::Power(sigma,7);
2551 fMu8 /= TMath::Power(sigma,8);
1a879e25 2552
f0754617 2553 // ----------------------------------------------------------------------
a123fe62 2554 // then calculate the higher moment errors
4e023902 2555 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2556 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2557 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
a123fe62 2558
2559 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2560
2561 // // assuming normal distribution (as done in ROOT)
2562 // meanError = sigma / TMath::Sqrt(normError);
2563 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2564 // skewnessError = TMath::Sqrt(6./(normError));
2565 // kurtosisError = TMath::Sqrt(24./(normError));
2566
2567 // use delta theorem paper (Luo - arXiv:1109.0593v1)
a9f20288 2568 Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2569 Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2570 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 2571 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2572 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2573 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
a123fe62 2574
4e023902 2575 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2576 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
a123fe62 2577
a9f20288 2578 if (TMath::Sqrt(normError) != 0){
2579 meanError = sigma / TMath::Sqrt(normError);
1a879e25 2580 sigmaError = TMath::Sqrt(Lambda11);
2581 skewnessError = TMath::Sqrt(Lambda22);
2582 kurtosisError = TMath::Sqrt(Lambda33);
2583
2584 success = kTRUE;
2585
a9f20288 2586 }
1a879e25 2587 else success = kFALSE;
2588
f2e8af26 2589 }
f2e8af26 2590 return success;
2591}
2592
2593
73609a48 2594//____________________________________________________________________//
2595Float_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) {
2596 //
2597 // calculates dphistar
2598 //
2599 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2600
2601 static const Double_t kPi = TMath::Pi();
2602
2603 // circularity
2604// if (dphistar > 2 * kPi)
2605// dphistar -= 2 * kPi;
2606// if (dphistar < -2 * kPi)
2607// dphistar += 2 * kPi;
2608
2609 if (dphistar > kPi)
2610 dphistar = kPi * 2 - dphistar;
2611 if (dphistar < -kPi)
2612 dphistar = -kPi * 2 - dphistar;
2613 if (dphistar > kPi) // might look funny but is needed
2614 dphistar = kPi * 2 - dphistar;
2615
2616 return dphistar;
2617}
2618
ea9867da 2619//____________________________________________________________________//
2620Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2621{
2622 // This method is a copy from AliUEHist::GetBinning
2623 // takes the binning from <configuration> identified by <tag>
2624 // configuration syntax example:
2625 // 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
2626 // phi: .....
2627 //
2628 // returns bin edges which have to be deleted by the caller
2629
2630 TString config(configuration);
2631 TObjArray* lines = config.Tokenize("\n");
2632 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2633 {
2634 TString line(lines->At(i)->GetName());
2635 if (line.BeginsWith(TString(tag) + ":"))
2636 {
2637 line.Remove(0, strlen(tag) + 1);
2638 line.ReplaceAll(" ", "");
2639 TObjArray* binning = line.Tokenize(",");
2640 Double_t* bins = new Double_t[binning->GetEntriesFast()];
2641 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2642 bins[j] = TString(binning->At(j)->GetName()).Atof();
2643
2644 nBins = binning->GetEntriesFast() - 1;
2922bbe3 2645
ea9867da 2646 delete binning;
2647 delete lines;
2648 return bins;
2649 }
2650 }
2651
2652 delete lines;
2653 AliFatal(Form("Tag %s not found in %s", tag, configuration));
2654 return 0;
2655}