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