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