]> 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
Memory leak fixed (Ruben)
[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());
399 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
400 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
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];
643 fHistConversionbefore->Fill(deta,dphi);
644
645 Float_t m0 = 0.510e-3;
646 Float_t tantheta1 = 1e10;
647
648 // phi in rad
c8ad9635 649 //Float_t phi1rad = firstPhi*TMath::DegToRad();
650 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
651 Float_t phi1rad = firstPhi;
652 Float_t phi2rad = secondPhi[j];
73609a48 653
654 if (firstEta < -1e-10 || firstEta > 1e-10)
655 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
656
657 Float_t tantheta2 = 1e10;
658 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
659 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
660
661 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
662 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
663
664 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
665
a9737921 666 if (masssqu < fInvMassCutConversion*fInvMassCutConversion){
73609a48 667 //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));
668 continue;
669 }
670 fHistConversionafter->Fill(deta,dphi);
671 }
672 }//conversion cut
5a861dc6 673
674 // momentum difference cut - suppress femtoscopic effects
675 if(fQCut){
676
6fa567bd 677 //Double_t ptMin = 0.1; //const for the time being (should be changeable later on)
41a6ae06 678 Double_t ptDifference = TMath::Abs( firstPt - secondPt[j]);
5a861dc6 679
41a6ae06 680 fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
6fa567bd 681 if(ptDifference < fDeltaPtMin) continue;
41a6ae06 682 fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
5a861dc6 683
684 }
685
35aff0f3 686 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
687 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
688 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
689 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
8795e993 690 else {
691 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
692 continue;
693 }
6acdbcb2 694 }//end of 2nd particle loop
695 }//end of 1st particle loop
0879e280 696}
697
0879e280 698//____________________________________________________________________//
9afe3098 699TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
700 Int_t iVariablePair,
9afe3098 701 Double_t psiMin,
6acdbcb2 702 Double_t psiMax,
20006629 703 Double_t vertexZMin,
704 Double_t vertexZMax,
6acdbcb2 705 Double_t ptTriggerMin,
706 Double_t ptTriggerMax,
707 Double_t ptAssociatedMin,
708 Double_t ptAssociatedMax) {
9afe3098 709 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 710 //(private members) of the AliBalancePsi class.
621329de 711 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
712 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 713
714 // security checks
b0d53e32 715 if(psiMin > psiMax-0.00001){
716 AliError("psiMax <= psiMin");
717 return NULL;
718 }
20006629 719 if(vertexZMin > vertexZMax-0.00001){
720 AliError("vertexZMax <= vertexZMin");
721 return NULL;
722 }
b0d53e32 723 if(ptTriggerMin > ptTriggerMax-0.00001){
724 AliError("ptTriggerMax <= ptTriggerMin");
725 return NULL;
726 }
727 if(ptAssociatedMin > ptAssociatedMax-0.00001){
728 AliError("ptAssociatedMax <= ptAssociatedMin");
729 return NULL;
730 }
622473b9 731
9afe3098 732 // Psi_2
622473b9 733 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
734 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
735 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
736 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
737 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
738 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
9afe3098 739
20006629 740 // Vz
741 if(fVertexBinning){
742 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
743 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
744 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
745 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
746 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
747 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
748 }
749
6acdbcb2 750 // pt trigger
751 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 752 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
753 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
754 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
755 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
756 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
757 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 758 }
759
760 // pt associated
761 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 762 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
763 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
764 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
765 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 766 }
767
a38fd7f3 768 //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));
769
9afe3098 770 // Project into the wanted space (1st: analysis step, 2nd: axis)
20006629 771 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); //
772 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); //
773 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); //
774 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); //
775 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); //
776 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); //
9afe3098 777
778 TH1D *gHistBalanceFunctionHistogram = 0x0;
779 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
780 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
781 gHistBalanceFunctionHistogram->Reset();
782
783 switch(iVariablePair) {
621329de 784 case 1:
c8ad9635 785 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
786 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 787 break;
621329de 788 case 2:
c8ad9635 789 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
790 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 791 break;
9afe3098 792 default:
793 break;
794 }
0879e280 795
0879e280 796 hTemp1->Sumw2();
797 hTemp2->Sumw2();
798 hTemp3->Sumw2();
799 hTemp4->Sumw2();
a38fd7f3 800 hTemp1->Add(hTemp3,-1.);
0879e280 801 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 802 hTemp2->Add(hTemp4,-1.);
0879e280 803 hTemp2->Scale(1./hTemp6->GetEntries());
804 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 805 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 806
807 //normalize to bin width
a9f20288 808 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
0879e280 809 }
810
0879e280 811 return gHistBalanceFunctionHistogram;
812}
a38fd7f3 813
f0c5040c 814//____________________________________________________________________//
815TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
816 Int_t iVariablePair,
817 Double_t psiMin,
818 Double_t psiMax,
20006629 819 Double_t vertexZMin,
820 Double_t vertexZMax,
f0c5040c 821 Double_t ptTriggerMin,
822 Double_t ptTriggerMax,
823 Double_t ptAssociatedMin,
824 Double_t ptAssociatedMax,
825 AliBalancePsi *bfMix) {
826 //Returns the BF histogram, extracted from the 6 AliTHn objects
827 //after dividing each correlation function by the Event Mixing one
828 //(private members) of the AliBalancePsi class.
829 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
830 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 831
832 // security checks
b0d53e32 833 if(psiMin > psiMax-0.00001){
834 AliError("psiMax <= psiMin");
835 return NULL;
836 }
20006629 837 if(vertexZMin > vertexZMax-0.00001){
838 AliError("vertexZMax <= vertexZMin");
839 return NULL;
840 }
b0d53e32 841 if(ptTriggerMin > ptTriggerMax-0.00001){
842 AliError("ptTriggerMax <= ptTriggerMin");
843 return NULL;
844 }
845 if(ptAssociatedMin > ptAssociatedMax-0.00001){
846 AliError("ptAssociatedMax <= ptAssociatedMin");
847 return NULL;
848 }
20006629 849
34616eda 850 // pt trigger (fixed for all vertices/psi/centralities)
f0c5040c 851 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 852 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
853 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
854 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
855 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
856 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
857 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 858 }
859
34616eda 860 // pt associated (fixed for all vertices/psi/centralities)
f0c5040c 861 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 862 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
863 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
864 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
865 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 866 }
867
868 // ============================================================================================
869 // the same for event mixing
870 AliTHn *fHistPMix = bfMix->GetHistNp();
871 AliTHn *fHistNMix = bfMix->GetHistNn();
872 AliTHn *fHistPNMix = bfMix->GetHistNpn();
873 AliTHn *fHistNPMix = bfMix->GetHistNnp();
874 AliTHn *fHistPPMix = bfMix->GetHistNpp();
875 AliTHn *fHistNNMix = bfMix->GetHistNnn();
876
34616eda 877 // pt trigger (fixed for all vertices/psi/centralities)
f0c5040c 878 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 879 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
880 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
881 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
882 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
883 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
884 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 885 }
886
34616eda 887 // pt associated (fixed for all vertices/psi/centralities)
f0c5040c 888 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 889 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
890 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
891 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
892 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 893 }
34616eda 894
f0c5040c 895 // ============================================================================================
34616eda 896 // ranges (in bins) for vertexZ and centrality (psi)
897 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
898 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
a4e2c68f 899 Int_t binVertexMin = 0;
900 Int_t binVertexMax = 0;
901 if(fVertexBinning){
902 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
903 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
904 }
34616eda 905 Double_t binPsiLowEdge = 0.;
a4e2c68f 906 Double_t binPsiUpEdge = 1000.;
34616eda 907 Double_t binVertexLowEdge = 0.;
a4e2c68f 908 Double_t binVertexUpEdge = 1000.;
34616eda 909
910 TH1D* h1 = NULL;
911 TH1D* h2 = NULL;
912 TH1D* h3 = NULL;
913 TH1D* h4 = NULL;
914
915 // loop over all bins
916 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
917 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
918
919 cout<<"In the balance function (1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
f0c5040c 920
34616eda 921 // determine the bin edges for this bin
922 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
923 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
a4e2c68f 924 if(fVertexBinning){
925 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
926 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
927 }
34616eda 928
929 // Psi_2
930 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
931 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
932 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
933 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
934 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
935 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
936
937 // Vz
938 if(fVertexBinning){
939 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
940 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
941 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
942 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
943 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
944 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
945 }
f0c5040c 946
34616eda 947 // ============================================================================================
948 // the same for event mixing
949
950 // Psi_2
951 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
952 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
953 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
954 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
955 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
956 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
957
958 // Vz
959 if(fVertexBinning){
960 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
961 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
962 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
963 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
964 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
965 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
966 }
967
968 // ============================================================================================
f0c5040c 969
34616eda 970 //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));
971
972 // Project into the wanted space (1st: analysis step, 2nd: axis)
973 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
974 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
975 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
976 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
977 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
978 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
979
980 // ============================================================================================
981 // the same for event mixing
982 TH1D* hTemp1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
983 TH1D* hTemp2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
984 TH1D* hTemp3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
985 TH1D* hTemp4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
986 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
987 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
988 // ============================================================================================
989
990 hTemp1->Sumw2();
991 hTemp2->Sumw2();
992 hTemp3->Sumw2();
993 hTemp4->Sumw2();
994 hTemp1Mix->Sumw2();
995 hTemp2Mix->Sumw2();
996 hTemp3Mix->Sumw2();
997 hTemp4Mix->Sumw2();
998
999 hTemp1->Scale(1./hTemp5->GetEntries());
1000 hTemp3->Scale(1./hTemp5->GetEntries());
1001 hTemp2->Scale(1./hTemp6->GetEntries());
1002 hTemp4->Scale(1./hTemp6->GetEntries());
271ceae3 1003
1004 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1005 // does not work here, so normalize also to trigger particles
1006 // --> careful: gives different integrals then as with full 2D method
34616eda 1007 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1008 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1009 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1010 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
271ceae3 1011
34616eda 1012 hTemp1->Divide(hTemp1Mix);
1013 hTemp2->Divide(hTemp2Mix);
1014 hTemp3->Divide(hTemp3Mix);
1015 hTemp4->Divide(hTemp4Mix);
1016
1017 // for the first: clone
1018 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1019 h1 = (TH1D*)hTemp1->Clone();
1020 h2 = (TH1D*)hTemp2->Clone();
1021 h3 = (TH1D*)hTemp3->Clone();
1022 h4 = (TH1D*)hTemp4->Clone();
1023 }
1024 else{ // otherwise: add for averaging
1025 h1->Add(hTemp1);
1026 h2->Add(hTemp2);
1027 h3->Add(hTemp3);
1028 h4->Add(hTemp4);
1029 }
1030
1031 }
1032 }
f0c5040c 1033
1034 TH1D *gHistBalanceFunctionHistogram = 0x0;
34616eda 1035 if((h1)&&(h2)&&(h3)&&(h4)) {
1036
1037 // average over number of bins nbinsVertex * nbinsPsi
1038 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1039 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1040 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1041 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1042
1043 gHistBalanceFunctionHistogram = (TH1D*)h1->Clone();
f0c5040c 1044 gHistBalanceFunctionHistogram->Reset();
1045
1046 switch(iVariablePair) {
1047 case 1:
1048 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1049 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1050 break;
1051 case 2:
1052 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1053 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1054 break;
1055 default:
1056 break;
1057 }
1058
34616eda 1059 h1->Add(h3,-1.);
1060 h2->Add(h4,-1.);
f0c5040c 1061
34616eda 1062 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
f0c5040c 1063 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1064
1065 //normalize to bin width
a9f20288 1066 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
f0c5040c 1067 }
1068
1069 return gHistBalanceFunctionHistogram;
1070}
1071
a38fd7f3 1072//____________________________________________________________________//
621329de 1073TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 1074 Double_t psiMax,
20006629 1075 Double_t vertexZMin,
1076 Double_t vertexZMax,
6acdbcb2 1077 Double_t ptTriggerMin,
1078 Double_t ptTriggerMax,
1079 Double_t ptAssociatedMin,
1080 Double_t ptAssociatedMax) {
621329de 1081 //Returns the BF histogram in Delta eta vs Delta phi,
1082 //extracted from the 6 AliTHn objects
1083 //(private members) of the AliBalancePsi class.
1084 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1085 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 1086
1087 // security checks
b0d53e32 1088 if(psiMin > psiMax-0.00001){
1089 AliError("psiMax <= psiMin");
1090 return NULL;
1091 }
20006629 1092 if(vertexZMin > vertexZMax-0.00001){
1093 AliError("vertexZMax <= vertexZMin");
1094 return NULL;
1095 }
b0d53e32 1096 if(ptTriggerMin > ptTriggerMax-0.00001){
1097 AliError("ptTriggerMax <= ptTriggerMin");
1098 return NULL;
1099 }
1100 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1101 AliError("ptAssociatedMax <= ptAssociatedMin");
1102 return NULL;
1103 }
622473b9 1104
621329de 1105 TString histName = "gHistBalanceFunctionHistogram2D";
1106
1107 // Psi_2
622473b9 1108 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1109 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1110 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1111 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1112 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1113 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1114
20006629 1115 // Vz
1116 if(fVertexBinning){
1117 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1118 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1119 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1120 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1121 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1122 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1123 }
1124
6acdbcb2 1125 // pt trigger
1126 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1127 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1128 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1129 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1130 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1131 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1132 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1133 }
1134
1135 // pt associated
1136 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1137 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1138 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1139 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1140 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1141 }
1142
1143 //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 1144
1145 // Project into the wanted space (1st: analysis step, 2nd: axis)
1146 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1147 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1148 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1149 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1150 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1151 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1152
1153 TH2D *gHistBalanceFunctionHistogram = 0x0;
1154 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1155 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1156 gHistBalanceFunctionHistogram->Reset();
c8ad9635 1157 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1158 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1159 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 1160
1161 hTemp1->Sumw2();
1162 hTemp2->Sumw2();
1163 hTemp3->Sumw2();
1164 hTemp4->Sumw2();
1165 hTemp1->Add(hTemp3,-1.);
1166 hTemp1->Scale(1./hTemp5->GetEntries());
1167 hTemp2->Add(hTemp4,-1.);
1168 hTemp2->Scale(1./hTemp6->GetEntries());
1169 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1170 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1171
1172 //normalize to bin width
1173 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
621329de 1174 }
1175
1176 return gHistBalanceFunctionHistogram;
1177}
1178
f0c5040c 1179//____________________________________________________________________//
1180TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1181 Double_t psiMax,
20006629 1182 Double_t vertexZMin,
1183 Double_t vertexZMax,
f0c5040c 1184 Double_t ptTriggerMin,
1185 Double_t ptTriggerMax,
1186 Double_t ptAssociatedMin,
1187 Double_t ptAssociatedMax,
1188 AliBalancePsi *bfMix) {
1189 //Returns the BF histogram in Delta eta vs Delta phi,
1190 //after dividing each correlation function by the Event Mixing one
1191 //extracted from the 6 AliTHn objects
1192 //(private members) of the AliBalancePsi class.
1193 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1194 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1195 TString histName = "gHistBalanceFunctionHistogram2D";
1196
1197 if(!bfMix){
1198 AliError("balance function object for event mixing not available");
1199 return NULL;
1200 }
1201
622473b9 1202 // security checks
b0d53e32 1203 if(psiMin > psiMax-0.00001){
1204 AliError("psiMax <= psiMin");
1205 return NULL;
1206 }
20006629 1207 if(vertexZMin > vertexZMax-0.00001){
1208 AliError("vertexZMax <= vertexZMin");
1209 return NULL;
1210 }
b0d53e32 1211 if(ptTriggerMin > ptTriggerMax-0.00001){
1212 AliError("ptTriggerMax <= ptTriggerMin");
1213 return NULL;
1214 }
1215 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1216 AliError("ptAssociatedMax <= ptAssociatedMin");
1217 return NULL;
1218 }
622473b9 1219
34616eda 1220 // pt trigger (fixed for all vertices/psi/centralities)
f0c5040c 1221 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1222 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1223 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1224 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1225 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1226 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1227 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1228 }
1229
34616eda 1230 // pt associated (fixed for all vertices/psi/centralities)
f0c5040c 1231 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1232 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1233 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1234 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1235 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1236 }
1237
f0c5040c 1238 // ============================================================================================
1239 // the same for event mixing
1240 AliTHn *fHistPMix = bfMix->GetHistNp();
1241 AliTHn *fHistNMix = bfMix->GetHistNn();
1242 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1243 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1244 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1245 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1246
34616eda 1247 // pt trigger (fixed for all vertices/psi/centralities)
f0c5040c 1248 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1249 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1250 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1251 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1252 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1253 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1254 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1255 }
1256
34616eda 1257 // pt associated (fixed for all vertices/psi/centralities)
f0c5040c 1258 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1259 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1260 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1261 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1262 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1263 }
34616eda 1264
f0c5040c 1265 // ============================================================================================
34616eda 1266 // ranges (in bins) for vertexZ and centrality (psi)
1267 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1268 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
a4e2c68f 1269 Int_t binVertexMin = 0;
1270 Int_t binVertexMax = 0;
1271 if(fVertexBinning){
1272 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1273 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1274 }
34616eda 1275 Double_t binPsiLowEdge = 0.;
1276 Double_t binPsiUpEdge = 0.;
1277 Double_t binVertexLowEdge = 0.;
1278 Double_t binVertexUpEdge = 0.;
1279
1280 TH2D* h1 = NULL;
1281 TH2D* h2 = NULL;
1282 TH2D* h3 = NULL;
1283 TH2D* h4 = NULL;
1284
1285 // loop over all bins
1286 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1287 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1288
1289 cout<<"In the balance function (2D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1290
1291 // determine the bin edges for this bin
1292 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1293 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
a4e2c68f 1294 if(fVertexBinning){
1295 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1296 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1297 }
34616eda 1298
1299 // Psi_2
1300 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1301 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1302 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1303 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1304 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1305 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1306
1307 // Vz
1308 if(fVertexBinning){
1309 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1310 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1311 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1312 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1313 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1314 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1315 }
f0c5040c 1316
1317
f0c5040c 1318
34616eda 1319 // ============================================================================================
1320 // the same for event mixing
1321
1322 // Psi_2
1323 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1324 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1325 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1326 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1327 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1328 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1329
1330 // Vz
1331 if(fVertexBinning){
1332 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1333 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1334 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1335 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1336 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1337 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1338 }
1339
1340 // ============================================================================================
1341
1342
1343 //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)));
1344
1345 // Project into the wanted space (1st: analysis step, 2nd: axis)
1346 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1347 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1348 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1349 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1350 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1351 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1352
1353 // ============================================================================================
1354 // the same for event mixing
1355 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1356 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1357 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1358 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
d631c24b 1359 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1360 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
34616eda 1361 // ============================================================================================
1362
1363 hTemp1->Sumw2();
1364 hTemp2->Sumw2();
1365 hTemp3->Sumw2();
1366 hTemp4->Sumw2();
1367 hTemp1Mix->Sumw2();
1368 hTemp2Mix->Sumw2();
1369 hTemp3Mix->Sumw2();
1370 hTemp4Mix->Sumw2();
1371
1372 hTemp1->Scale(1./hTemp5->GetEntries());
1373 hTemp3->Scale(1./hTemp5->GetEntries());
1374 hTemp2->Scale(1./hTemp6->GetEntries());
1375 hTemp4->Scale(1./hTemp6->GetEntries());
271ceae3 1376
1377 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1378 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1379 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1380 hTemp1Mix->Scale(1./mixedNorm1);
1381 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1382 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1383 hTemp2Mix->Scale(1./mixedNorm2);
1384 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1385 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1386 hTemp3Mix->Scale(1./mixedNorm3);
1387 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1388 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1389 hTemp4Mix->Scale(1./mixedNorm4);
34616eda 1390
1391 hTemp1->Divide(hTemp1Mix);
1392 hTemp2->Divide(hTemp2Mix);
1393 hTemp3->Divide(hTemp3Mix);
1394 hTemp4->Divide(hTemp4Mix);
1395
1396 // for the first: clone
1397 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1398 h1 = (TH2D*)hTemp1->Clone();
1399 h2 = (TH2D*)hTemp2->Clone();
1400 h3 = (TH2D*)hTemp3->Clone();
1401 h4 = (TH2D*)hTemp4->Clone();
1402 }
1403 else{ // otherwise: add for averaging
1404 h1->Add(hTemp1);
1405 h2->Add(hTemp2);
1406 h3->Add(hTemp3);
1407 h4->Add(hTemp4);
1408 }
1409 }
1410 }
1411
1412 TH2D *gHistBalanceFunctionHistogram = 0x0;
1413 if((h1)&&(h2)&&(h3)&&(h4)) {
f0c5040c 1414
34616eda 1415 // average over number of bins nbinsVertex * nbinsPsi
1416 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1417 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1418 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1419 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
f0c5040c 1420
34616eda 1421 gHistBalanceFunctionHistogram = (TH2D*)h1->Clone();
f0c5040c 1422 gHistBalanceFunctionHistogram->Reset();
1423 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1424 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1425 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
f0c5040c 1426
34616eda 1427 h1->Add(h3,-1.);
1428 h2->Add(h4,-1.);
1429
1430 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
f0c5040c 1431 gHistBalanceFunctionHistogram->Scale(0.5);
34616eda 1432
faa03677 1433 //normalize to bin width
1434 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 1435 }
34616eda 1436
f0c5040c 1437 return gHistBalanceFunctionHistogram;
1438}
1439
f2e8af26 1440//____________________________________________________________________//
1441TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1442 Double_t psiMin,
1443 Double_t psiMax,
20006629 1444 Double_t vertexZMin,
1445 Double_t vertexZMax,
f2e8af26 1446 Double_t ptTriggerMin,
1447 Double_t ptTriggerMax,
1448 Double_t ptAssociatedMin,
1449 Double_t ptAssociatedMax,
1450 AliBalancePsi *bfMix) {
1451 //Returns the BF histogram in Delta eta OR Delta phi,
1452 //after dividing each correlation function by the Event Mixing one
1453 // (But the division is done here in 2D, this was basically done to check the Integral)
1454 //extracted from the 6 AliTHn objects
1455 //(private members) of the AliBalancePsi class.
1456 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1457 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1458 TString histName = "gHistBalanceFunctionHistogram1D";
1459
1460 if(!bfMix){
1461 AliError("balance function object for event mixing not available");
1462 return NULL;
1463 }
1464
622473b9 1465 // security checks
b0d53e32 1466 if(psiMin > psiMax-0.00001){
1467 AliError("psiMax <= psiMin");
1468 return NULL;
1469 }
20006629 1470 if(vertexZMin > vertexZMax-0.00001){
1471 AliError("vertexZMax <= vertexZMin");
1472 return NULL;
1473 }
b0d53e32 1474 if(ptTriggerMin > ptTriggerMax-0.00001){
1475 AliError("ptTriggerMax <= ptTriggerMin");
1476 return NULL;
1477 }
1478 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1479 AliError("ptAssociatedMax <= ptAssociatedMin");
1480 return NULL;
1481 }
622473b9 1482
34616eda 1483 // pt trigger (fixed for all vertices/psi/centralities)
f2e8af26 1484 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1485 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1486 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1487 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1488 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1489 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1490 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1491 }
1492
34616eda 1493 // pt associated (fixed for all vertices/psi/centralities)
f2e8af26 1494 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1495 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1496 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1497 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1498 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1499 }
1500
f2e8af26 1501 // ============================================================================================
1502 // the same for event mixing
1503 AliTHn *fHistPMix = bfMix->GetHistNp();
1504 AliTHn *fHistNMix = bfMix->GetHistNn();
1505 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1506 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1507 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1508 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1509
34616eda 1510 // pt trigger (fixed for all vertices/psi/centralities)
f2e8af26 1511 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1512 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1513 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1514 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1515 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1516 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1517 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1518 }
1519
34616eda 1520 // pt associated (fixed for all vertices/psi/centralities)
f2e8af26 1521 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1522 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1523 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1524 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1525 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1526 }
f2e8af26 1527
34616eda 1528 // ============================================================================================
1529 // ranges (in bins) for vertexZ and centrality (psi)
1530 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1531 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
a4e2c68f 1532 Int_t binVertexMin = 0;
1533 Int_t binVertexMax = 0;
1534 if(fVertexBinning){
1535 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1536 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1537 }
34616eda 1538 Double_t binPsiLowEdge = 0.;
1539 Double_t binPsiUpEdge = 0.;
1540 Double_t binVertexLowEdge = 0.;
1541 Double_t binVertexUpEdge = 0.;
1542
1543 TH2D* h1 = NULL;
1544 TH2D* h2 = NULL;
1545 TH2D* h3 = NULL;
1546 TH2D* h4 = NULL;
1547
1548 // loop over all bins
1549 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1550 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1551
1552 cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1553
1554 // determine the bin edges for this bin
1555 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1556 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
a4e2c68f 1557 if(fVertexBinning){
1558 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1559 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1560 }
34616eda 1561
1562 // Psi_2
1563 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1564 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1565 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1566 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1567 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1568 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1569
1570 // Vz
1571 if(fVertexBinning){
1572 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1573 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1574 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1575 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1576 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1577 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1578 }
f2e8af26 1579
34616eda 1580 // ============================================================================================
1581 // the same for event mixing
1582
1583 // Psi_2
1584 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1585 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1586 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1587 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1588 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1589 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1590
1591 // Vz
1592 if(fVertexBinning){
1593 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1594 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1595 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1596 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1597 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1598 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1599 }
1600 // ============================================================================================
f2e8af26 1601
34616eda 1602 //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)));
1603
1604 // Project into the wanted space (1st: analysis step, 2nd: axis)
1605 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1606 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1607 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1608 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1609 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1610 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1611
1612 // ============================================================================================
1613 // the same for event mixing
1614 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1615 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1616 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1617 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
d631c24b 1618 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1619 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
34616eda 1620 // ============================================================================================
1621
1622 hTemp1->Sumw2();
1623 hTemp2->Sumw2();
1624 hTemp3->Sumw2();
1625 hTemp4->Sumw2();
1626 hTemp1Mix->Sumw2();
1627 hTemp2Mix->Sumw2();
1628 hTemp3Mix->Sumw2();
1629 hTemp4Mix->Sumw2();
1630
1631 hTemp1->Scale(1./hTemp5->GetEntries());
1632 hTemp3->Scale(1./hTemp5->GetEntries());
1633 hTemp2->Scale(1./hTemp6->GetEntries());
1634 hTemp4->Scale(1./hTemp6->GetEntries());
271ceae3 1635
1636 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1637 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1638 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1639 hTemp1Mix->Scale(1./mixedNorm1);
1640 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1641 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1642 hTemp2Mix->Scale(1./mixedNorm2);
1643 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1644 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1645 hTemp3Mix->Scale(1./mixedNorm3);
1646 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1647 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1648 hTemp4Mix->Scale(1./mixedNorm4);
1649
34616eda 1650 hTemp1->Divide(hTemp1Mix);
1651 hTemp2->Divide(hTemp2Mix);
1652 hTemp3->Divide(hTemp3Mix);
1653 hTemp4->Divide(hTemp4Mix);
1654
1655 // for the first: clone
1656 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1657 h1 = (TH2D*)hTemp1->Clone();
1658 h2 = (TH2D*)hTemp2->Clone();
1659 h3 = (TH2D*)hTemp3->Clone();
1660 h4 = (TH2D*)hTemp4->Clone();
1661 }
1662 else{ // otherwise: add for averaging
1663 h1->Add(hTemp1);
1664 h2->Add(hTemp2);
1665 h3->Add(hTemp3);
1666 h4->Add(hTemp4);
1667 }
f2e8af26 1668
34616eda 1669 }
1670 }
f2e8af26 1671
1672 TH1D *gHistBalanceFunctionHistogram = 0x0;
34616eda 1673 if((h1)&&(h2)&&(h3)&&(h4)) {
f2e8af26 1674
34616eda 1675 // average over number of bins nbinsVertex * nbinsPsi
1676 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1677 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1678 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1679 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
f2e8af26 1680
1681 // now only project on one axis
1682 TH1D *h1DTemp1 = NULL;
1683 TH1D *h1DTemp2 = NULL;
1684 TH1D *h1DTemp3 = NULL;
1685 TH1D *h1DTemp4 = NULL;
1686 if(!bPhi){
34616eda 1687 h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName()));
1688 h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName()));
1689 h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName()));
1690 h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName()));
f2e8af26 1691 }
1692 else{
34616eda 1693 h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName()));
1694 h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName()));
1695 h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName()));
1696 h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName()));
f2e8af26 1697 }
1698
1699 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1700 gHistBalanceFunctionHistogram->Reset();
1701 if(!bPhi){
1702 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1703 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1704 }
1705 else{
1706 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1707 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1708 }
1709
1710 h1DTemp1->Add(h1DTemp3,-1.);
1711 h1DTemp2->Add(h1DTemp4,-1.);
1712
1713 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1714 gHistBalanceFunctionHistogram->Scale(0.5);
a9f20288 1715
1716 //normalize to bin width
1717 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
f2e8af26 1718 }
1719
1720 return gHistBalanceFunctionHistogram;
1721}
1722
f0c5040c 1723
34616eda 1724//____________________________________________________________________//
1725TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1726 Double_t psiMin,
1727 Double_t psiMax,
1728 Double_t vertexZMin,
1729 Double_t vertexZMax,
1730 Double_t ptTriggerMin,
1731 Double_t ptTriggerMax,
1732 Double_t ptAssociatedMin,
1733 Double_t ptAssociatedMax,
1734 AliBalancePsi *bMixed) {
1735
1736 // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1737 // does the division by event mixing inside,
1738 // and averaging over several vertexZ and centrality bins
1739
1740 // security checks
1741 if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1742 AliError("Only types allowed: PN,NP,PP,NN,ALL");
1743 return NULL;
1744 }
1745 if(!bMixed){
1746 AliError("No Event Mixing AliTHn");
1747 return NULL;
1748 }
1749
1750 TH2D *gHist = NULL;
1751 TH2D *fSame = NULL;
1752 TH2D *fMixed = NULL;
1753
1754 // ranges (in bins) for vertexZ and centrality (psi)
1755 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1756 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
a4e2c68f 1757 Int_t binVertexMin = 0;
1758 Int_t binVertexMax = 0;
1759 if(fVertexBinning){
1760 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1761 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1762 }
34616eda 1763 Double_t binPsiLowEdge = 0.;
a4e2c68f 1764 Double_t binPsiUpEdge = 1000.;
34616eda 1765 Double_t binVertexLowEdge = 0.;
a4e2c68f 1766 Double_t binVertexUpEdge = 1000.;
34616eda 1767
1768 // loop over all bins
1769 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1770 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1771
1772 cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1773
1774 // determine the bin edges for this bin
1775 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1776 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
a4e2c68f 1777 if(fVertexBinning){
1778 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1779 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1780 }
1781
34616eda 1782 // get the 2D histograms for the correct type
1783 if(type=="PN"){
1784 fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1785 fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1786 }
1787 else if(type=="NP"){
1788 fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1789 fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1790 }
1791 else if(type=="PP"){
1792 fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
7071dbe3 1793 fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
34616eda 1794 }
1795 else if(type=="NN"){
1796 fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1797 fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1798 }
1799 else if(type=="ALL"){
1800 fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1801 fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1802 }
1803
1fdb6a72 1804 if(fSame && fMixed){
1805 // then get the correlation function (divide fSame/fmixed)
1806 fSame->Divide(fMixed);
1807
7071dbe3 1808 // NEW averaging:
1809 // average over number of triggers in each sub-bin
a9737921 1810 Double_t NTrigSubBin = 0;
1811 if(type=="PN" || type=="PP")
1812 NTrigSubBin = (Double_t)(fHistP->Project(0,1)->GetEntries());
1813 else if(type=="NP" || type=="NN")
1814 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->GetEntries());
7071dbe3 1815 fSame->Scale(NTrigSubBin);
1816
1fdb6a72 1817 // for the first: clone
84d14827 1818 if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
1fdb6a72 1819 gHist = (TH2D*)fSame->Clone();
1820 }
1821 else{ // otherwise: add for averaging
1822 gHist->Add(fSame);
1823 }
34616eda 1824 }
1825 }
1826 }
1827
1fdb6a72 1828 if(gHist){
7071dbe3 1829
1830 // OLD averaging:
1fdb6a72 1831 // average over number of bins nbinsVertex * nbinsPsi
7071dbe3 1832 // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1833
1834 // NEW averaging:
1835 // average over number of triggers in each sub-bin
1836 // first set to full range and then obtain number of all triggers
a9737921 1837 Double_t NTrigAll = 0;
1838 if(type=="PN" || type=="PP"){
1839 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1840 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1841 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1842 NTrigAll = (Double_t)(fHistP->Project(0,1)->GetEntries());
1843 }
1844 else if(type=="NP" || type=="NN"){
1845 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1846 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1847 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1848 NTrigAll = (Double_t)(fHistN->Project(0,1)->GetEntries());
1849 }
7071dbe3 1850 gHist->Scale(1./NTrigAll);
1851
1fdb6a72 1852 }
34616eda 1853
1854 return gHist;
1855}
1856
1857
621329de 1858//____________________________________________________________________//
1859TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1860 Double_t psiMax,
20006629 1861 Double_t vertexZMin,
1862 Double_t vertexZMax,
6acdbcb2 1863 Double_t ptTriggerMin,
1864 Double_t ptTriggerMax,
1865 Double_t ptAssociatedMin,
47e040b5 1866 Double_t ptAssociatedMax) {
a38fd7f3 1867 //Returns the 2D correlation function for (+-) pairs
622473b9 1868
1869 // security checks
b0d53e32 1870 if(psiMin > psiMax-0.00001){
1871 AliError("psiMax <= psiMin");
1872 return NULL;
1873 }
20006629 1874 if(vertexZMin > vertexZMax-0.00001){
1875 AliError("vertexZMax <= vertexZMin");
1876 return NULL;
1877 }
b0d53e32 1878 if(ptTriggerMin > ptTriggerMax-0.00001){
1879 AliError("ptTriggerMax <= ptTriggerMin");
1880 return NULL;
1881 }
1882 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1883 AliError("ptAssociatedMax <= ptAssociatedMin");
1884 return NULL;
1885 }
622473b9 1886
621329de 1887 // Psi_2: axis 0
622473b9 1888 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1889 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1890
20006629 1891 // Vz
1892 if(fVertexBinning){
1893 //Printf("Cutting on Vz...");
1894 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1895 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1896 }
1897
6acdbcb2 1898 // pt trigger
1899 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1900 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1901 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1902 }
621329de 1903
6acdbcb2 1904 // pt associated
1905 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1906 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1907
1908 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1909 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1910
1911 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1912 //TCanvas *c1 = new TCanvas("c1","");
1913 //c1->cd();
1914 //if(!gHistTest){
1915 //AliError("Projection of fHistP = NULL");
1916 //return gHistTest;
1917 //}
1918 //else{
1919 //gHistTest->DrawCopy("colz");
1920 //}
1921
1922 //0:step, 1: Delta eta, 2: Delta phi
621329de 1923 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1924 if(!gHist){
1925 AliError("Projection of fHistPN = NULL");
1926 return gHist;
1927 }
1928
6acdbcb2 1929 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1930 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1931 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1932
6acdbcb2 1933 //TCanvas *c2 = new TCanvas("c2","");
1934 //c2->cd();
1935 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1936
47e040b5 1937 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1938 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1939
1940 //normalize to bin width
1941 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1942
1943 return gHist;
1944}
1945
1946//____________________________________________________________________//
621329de 1947TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1948 Double_t psiMax,
20006629 1949 Double_t vertexZMin,
1950 Double_t vertexZMax,
6acdbcb2 1951 Double_t ptTriggerMin,
1952 Double_t ptTriggerMax,
1953 Double_t ptAssociatedMin,
47e040b5 1954 Double_t ptAssociatedMax) {
a38fd7f3 1955 //Returns the 2D correlation function for (+-) pairs
622473b9 1956
1957 // security checks
b0d53e32 1958 if(psiMin > psiMax-0.00001){
1959 AliError("psiMax <= psiMin");
1960 return NULL;
1961 }
20006629 1962 if(vertexZMin > vertexZMax-0.00001){
1963 AliError("vertexZMax <= vertexZMin");
1964 return NULL;
1965 }
b0d53e32 1966 if(ptTriggerMin > ptTriggerMax-0.00001){
1967 AliError("ptTriggerMax <= ptTriggerMin");
1968 return NULL;
1969 }
1970 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1971 AliError("ptAssociatedMax <= ptAssociatedMin");
1972 return NULL;
1973 }
622473b9 1974
621329de 1975 // Psi_2: axis 0
622473b9 1976 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1977 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 1978
1979 // Vz
1980 if(fVertexBinning){
1981 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1982 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1983 }
a38fd7f3 1984
6acdbcb2 1985 // pt trigger
1986 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1987 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1988 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1989 }
1990
1991 // pt associated
1992 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1993 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1994
1995 //0:step, 1: Delta eta, 2: Delta phi
621329de 1996 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1997 if(!gHist){
1998 AliError("Projection of fHistPN = NULL");
1999 return gHist;
2000 }
2001
a38fd7f3 2002 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2003 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 2004 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 2005 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 2006
2007 //normalize to bin width
2008 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2009
2010 return gHist;
2011}
2012
2013//____________________________________________________________________//
621329de 2014TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 2015 Double_t psiMax,
20006629 2016 Double_t vertexZMin,
2017 Double_t vertexZMax,
6acdbcb2 2018 Double_t ptTriggerMin,
2019 Double_t ptTriggerMax,
2020 Double_t ptAssociatedMin,
47e040b5 2021 Double_t ptAssociatedMax) {
a38fd7f3 2022 //Returns the 2D correlation function for (+-) pairs
622473b9 2023
2024 // security checks
b0d53e32 2025 if(psiMin > psiMax-0.00001){
2026 AliError("psiMax <= psiMin");
2027 return NULL;
2028 }
20006629 2029 if(vertexZMin > vertexZMax-0.00001){
2030 AliError("vertexZMax <= vertexZMin");
2031 return NULL;
2032 }
b0d53e32 2033 if(ptTriggerMin > ptTriggerMax-0.00001){
2034 AliError("ptTriggerMax <= ptTriggerMin");
2035 return NULL;
2036 }
2037 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2038 AliError("ptAssociatedMax <= ptAssociatedMin");
2039 return NULL;
2040 }
622473b9 2041
621329de 2042 // Psi_2: axis 0
622473b9 2043 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2044 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 2045
20006629 2046 // Vz
2047 if(fVertexBinning){
2048 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2049 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2050 }
2051
6acdbcb2 2052 // pt trigger
2053 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2054 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2055 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2056 }
2057
2058 // pt associated
2059 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2060 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 2061
6acdbcb2 2062 //0:step, 1: Delta eta, 2: Delta phi
621329de 2063 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 2064 if(!gHist){
2065 AliError("Projection of fHistPN = NULL");
2066 return gHist;
2067 }
2068
a38fd7f3 2069 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2070 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 2071 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 2072 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 2073
2074 //normalize to bin width
2075 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2076
2077 return gHist;
2078}
2079
2080//____________________________________________________________________//
621329de 2081TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 2082 Double_t psiMax,
20006629 2083 Double_t vertexZMin,
2084 Double_t vertexZMax,
6acdbcb2 2085 Double_t ptTriggerMin,
2086 Double_t ptTriggerMax,
2087 Double_t ptAssociatedMin,
47e040b5 2088 Double_t ptAssociatedMax) {
a38fd7f3 2089 //Returns the 2D correlation function for (+-) pairs
622473b9 2090
2091 // security checks
2092 if(psiMin > psiMax-0.00001){
2093 AliError("psiMax <= psiMin");
2094 return NULL;
2095 }
20006629 2096 if(vertexZMin > vertexZMax-0.00001){
2097 AliError("vertexZMax <= vertexZMin");
2098 return NULL;
2099 }
622473b9 2100 if(ptTriggerMin > ptTriggerMax-0.00001){
2101 AliError("ptTriggerMax <= ptTriggerMin");
2102 return NULL;
2103 }
2104 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2105 AliError("ptAssociatedMax <= ptAssociatedMin");
2106 return NULL;
2107 }
2108
621329de 2109 // Psi_2: axis 0
622473b9 2110 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2111 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 2112
20006629 2113 // Vz
2114 if(fVertexBinning){
2115 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2116 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2117 }
2118
6acdbcb2 2119 // pt trigger
2120 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2121 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2122 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2123 }
2124
2125 // pt associated
2126 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2127 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 2128
6acdbcb2 2129 //0:step, 1: Delta eta, 2: Delta phi
621329de 2130 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 2131 if(!gHist){
2132 AliError("Projection of fHistPN = NULL");
2133 return gHist;
2134 }
2135
a38fd7f3 2136 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2137 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 2138 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 2139 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 2140
2141 //normalize to bin width
2142 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2143
2144 return gHist;
2145}
621329de 2146
7b610f27 2147//____________________________________________________________________//
2148TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2149 Double_t psiMax,
20006629 2150 Double_t vertexZMin,
2151 Double_t vertexZMax,
7b610f27 2152 Double_t ptTriggerMin,
2153 Double_t ptTriggerMax,
2154 Double_t ptAssociatedMin,
47e040b5 2155 Double_t ptAssociatedMax) {
7b610f27 2156 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 2157
2158 // security checks
b0d53e32 2159 if(psiMin > psiMax-0.00001){
2160 AliError("psiMax <= psiMin");
2161 return NULL;
2162 }
20006629 2163 if(vertexZMin > vertexZMax-0.00001){
2164 AliError("vertexZMax <= vertexZMin");
2165 return NULL;
2166 }
b0d53e32 2167 if(ptTriggerMin > ptTriggerMax-0.00001){
2168 AliError("ptTriggerMax <= ptTriggerMin");
2169 return NULL;
2170 }
2171 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2172 AliError("ptAssociatedMax <= ptAssociatedMin");
2173 return NULL;
2174 }
622473b9 2175
7b610f27 2176 // Psi_2: axis 0
622473b9 2177 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2178 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2179 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2180 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2181 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2182 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 2183
2184 // Vz
2185 if(fVertexBinning){
2186 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2187 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2188 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2189 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2190 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2191 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2192 }
7b610f27 2193
2194 // pt trigger
2195 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2196 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2197 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2198 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2199 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2200 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2201 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 2202 }
2203
2204 // pt associated
2205 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2206 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2207 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2208 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2209 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 2210
2211 //0:step, 1: Delta eta, 2: Delta phi
2212 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2213 if(!gHistNN){
2214 AliError("Projection of fHistNN = NULL");
2215 return gHistNN;
2216 }
2217 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2218 if(!gHistPP){
2219 AliError("Projection of fHistPP = NULL");
2220 return gHistPP;
2221 }
2222 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2223 if(!gHistNP){
2224 AliError("Projection of fHistNP = NULL");
2225 return gHistNP;
2226 }
2227 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2228 if(!gHistPN){
2229 AliError("Projection of fHistPN = NULL");
2230 return gHistPN;
2231 }
2232
2233 // sum all 2 particle histograms
2234 gHistNN->Add(gHistPP);
2235 gHistNN->Add(gHistNP);
2236 gHistNN->Add(gHistPN);
2237
47e040b5 2238 // divide by sum of + and - triggers
2239 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 2240 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 2241
2242 //normalize to bin width
2243 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 2244
2245 return gHistNN;
2246}
2247
f2e8af26 2248//____________________________________________________________________//
2249
1822002d 2250Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
f2e8af26 2251 Double_t &mean, Double_t &meanError,
2252 Double_t &sigma, Double_t &sigmaError,
2253 Double_t &skewness, Double_t &skewnessError,
2254 Double_t &kurtosis, Double_t &kurtosisError) {
2255 //
2256 // helper method to calculate the moments and errors of a TH1D anlytically
4e023902 2257 // fVariable = 1 for Delta eta (moments in full range)
2258 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
f2e8af26 2259 //
2260
2261 Bool_t success = kFALSE;
2262 mean = -1.;
2263 meanError = -1.;
2264 sigma = -1.;
2265 sigmaError = -1.;
2266 skewness = -1.;
2267 skewnessError = -1.;
2268 kurtosis = -1.;
2269 kurtosisError = -1.;
2270
2271 if(gHist){
2272
2273 // ----------------------------------------------------------------------
2274 // basic parameters of histogram
2275
2276 Int_t fNumberOfBins = gHist->GetNbinsX();
2277 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
2278
1a879e25 2279
2280 // ----------------------------------------------------------------------
2281 // ZYAM (for partially negative distributions)
2282 // --> we subtract always the minimum value
7071dbe3 2283 Double_t zeroYield = 0.;
2284 Double_t zeroYieldCur = -FLT_MAX;
2285 if(kUseZYAM){
4184be06 2286 for(Int_t iMin = 0; iMin<2; iMin++){
7071dbe3 2287 zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2288 zeroYield += zeroYieldCur;
2289 }
4184be06 2290 zeroYield /= 2.;
7071dbe3 2291 //zeroYield = gHist->GetMinimum();
2292 }
f2e8af26 2293 // ----------------------------------------------------------------------
2294 // first calculate the mean
2295
2296 Double_t fWeightedAverage = 0.;
f0754617 2297 Double_t fNormalization = 0.;
f2e8af26 2298
2299 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2300
4e023902 2301 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2302 if(fVariable == 2 &&
2303 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2304 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2305 continue;
2306 }
2307
1a879e25 2308 fWeightedAverage += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2309 fNormalization += (gHist->GetBinContent(i)-zeroYield);
f2e8af26 2310 }
2311
2312 mean = fWeightedAverage / fNormalization;
2313
f2e8af26 2314 // ----------------------------------------------------------------------
2315 // then calculate the higher moments
2316
a123fe62 2317 Double_t fMu = 0.;
2318 Double_t fMu2 = 0.;
2319 Double_t fMu3 = 0.;
2320 Double_t fMu4 = 0.;
2321 Double_t fMu5 = 0.;
2322 Double_t fMu6 = 0.;
2323 Double_t fMu7 = 0.;
2324 Double_t fMu8 = 0.;
f2e8af26 2325
2326 for(Int_t i = 1; i <= fNumberOfBins; i++) {
4e023902 2327
2328 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2329 if(fVariable == 2 &&
2330 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2331 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2332 continue;
2333 }
2334
1a879e25 2335 fMu += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2336 fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2337 fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2338 fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2339 fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2340 fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2341 fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2342 fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
f2e8af26 2343 }
e47eeabd 2344
2345 // normalize to bin entries!
2346 fMu /= fNormalization;
2347 fMu2 /= fNormalization;
2348 fMu3 /= fNormalization;
2349 fMu4 /= fNormalization;
2350 fMu5 /= fNormalization;
2351 fMu6 /= fNormalization;
2352 fMu7 /= fNormalization;
2353 fMu8 /= fNormalization;
2354
2355 sigma = TMath::Sqrt(fMu2);
2356 skewness = fMu3 / TMath::Power(sigma,3);
2357 kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
f2e8af26 2358
efef044a 2359 // normalize with sigma^r
2360 fMu /= TMath::Power(sigma,1);
2361 fMu2 /= TMath::Power(sigma,2);
2362 fMu3 /= TMath::Power(sigma,3);
2363 fMu4 /= TMath::Power(sigma,4);
2364 fMu5 /= TMath::Power(sigma,5);
2365 fMu6 /= TMath::Power(sigma,6);
2366 fMu7 /= TMath::Power(sigma,7);
2367 fMu8 /= TMath::Power(sigma,8);
1a879e25 2368
f0754617 2369 // ----------------------------------------------------------------------
a123fe62 2370 // then calculate the higher moment errors
4e023902 2371 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2372 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2373 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
a123fe62 2374
2375 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2376
2377 // // assuming normal distribution (as done in ROOT)
2378 // meanError = sigma / TMath::Sqrt(normError);
2379 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2380 // skewnessError = TMath::Sqrt(6./(normError));
2381 // kurtosisError = TMath::Sqrt(24./(normError));
2382
2383 // use delta theorem paper (Luo - arXiv:1109.0593v1)
a9f20288 2384 Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2385 Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2386 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 2387 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2388 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2389 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
a123fe62 2390
4e023902 2391 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2392 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
a123fe62 2393
a9f20288 2394 if (TMath::Sqrt(normError) != 0){
2395 meanError = sigma / TMath::Sqrt(normError);
1a879e25 2396 sigmaError = TMath::Sqrt(Lambda11);
2397 skewnessError = TMath::Sqrt(Lambda22);
2398 kurtosisError = TMath::Sqrt(Lambda33);
2399
2400 success = kTRUE;
2401
a9f20288 2402 }
1a879e25 2403 else success = kFALSE;
2404
f2e8af26 2405 }
f2e8af26 2406 return success;
2407}
2408
2409
73609a48 2410//____________________________________________________________________//
2411Float_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) {
2412 //
2413 // calculates dphistar
2414 //
2415 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2416
2417 static const Double_t kPi = TMath::Pi();
2418
2419 // circularity
2420// if (dphistar > 2 * kPi)
2421// dphistar -= 2 * kPi;
2422// if (dphistar < -2 * kPi)
2423// dphistar += 2 * kPi;
2424
2425 if (dphistar > kPi)
2426 dphistar = kPi * 2 - dphistar;
2427 if (dphistar < -kPi)
2428 dphistar = -kPi * 2 - dphistar;
2429 if (dphistar > kPi) // might look funny but is needed
2430 dphistar = kPi * 2 - dphistar;
2431
2432 return dphistar;
2433}
2434
2922bbe3 2435
2436