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