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