]> 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
from Ante Bilandzic:
[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 ------------------------------------
a9f20288 180 const Int_t kMultBins = 10;
f0fb4ac1 181 //A first rough attempt at four bins
a9f20288 182 Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80,100,100000};
f0fb4ac1 183 //----------------------------------------------------------
184
185 //--- Centrality Bins --------------------------------------
683cbfb5 186 const Int_t kNCentralityBins = 9;
a9f20288 187 const Int_t kNCentralityBinsVertex = 15;
683cbfb5 188 Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
a9f20288 189 Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,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;
a9f20288 269 const Int_t kNPtBinsVertex = 6;
683cbfb5 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.};
a9f20288 271 Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0,15.0};
683cbfb5 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
a9f20288 799 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->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
a9f20288 1047 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->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);
a9f20288 1664
1665 //normalize to bin width
1666 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
f2e8af26 1667 }
1668
1669 return gHistBalanceFunctionHistogram;
1670}
1671
f0c5040c 1672
34616eda 1673//____________________________________________________________________//
1674TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1675 Double_t psiMin,
1676 Double_t psiMax,
1677 Double_t vertexZMin,
1678 Double_t vertexZMax,
1679 Double_t ptTriggerMin,
1680 Double_t ptTriggerMax,
1681 Double_t ptAssociatedMin,
1682 Double_t ptAssociatedMax,
1683 AliBalancePsi *bMixed) {
1684
1685 // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1686 // does the division by event mixing inside,
1687 // and averaging over several vertexZ and centrality bins
1688
1689 // security checks
1690 if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1691 AliError("Only types allowed: PN,NP,PP,NN,ALL");
1692 return NULL;
1693 }
1694 if(!bMixed){
1695 AliError("No Event Mixing AliTHn");
1696 return NULL;
1697 }
1698
1699 TH2D *gHist = NULL;
1700 TH2D *fSame = NULL;
1701 TH2D *fMixed = NULL;
1702
1703 // ranges (in bins) for vertexZ and centrality (psi)
1704 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1705 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1706 Int_t binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1707 Int_t binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1708 Double_t binPsiLowEdge = 0.;
1709 Double_t binPsiUpEdge = 0.;
1710 Double_t binVertexLowEdge = 0.;
1711 Double_t binVertexUpEdge = 0.;
1712
1713 // loop over all bins
1714 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1715 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1716
1717 cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1718
1719 // determine the bin edges for this bin
1720 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1721 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1722 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1723 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1724
1725 // get the 2D histograms for the correct type
1726 if(type=="PN"){
1727 fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1728 fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1729 }
1730 else if(type=="NP"){
1731 fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1732 fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1733 }
1734 else if(type=="PP"){
1735 fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1736 fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1737 }
1738 else if(type=="NN"){
1739 fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1740 fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1741 }
1742 else if(type=="ALL"){
1743 fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1744 fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1745 }
1746
1fdb6a72 1747 if(fSame && fMixed){
1748 // then get the correlation function (divide fSame/fmixed)
1749 fSame->Divide(fMixed);
1750
1751 // for the first: clone
84d14827 1752 if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
1fdb6a72 1753 gHist = (TH2D*)fSame->Clone();
1754 }
1755 else{ // otherwise: add for averaging
1756 gHist->Add(fSame);
1757 }
34616eda 1758 }
1759 }
1760 }
1761
1fdb6a72 1762 if(gHist){
1763 // average over number of bins nbinsVertex * nbinsPsi
1764 gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1765 }
34616eda 1766
1767 return gHist;
1768}
1769
1770
621329de 1771//____________________________________________________________________//
1772TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1773 Double_t psiMax,
20006629 1774 Double_t vertexZMin,
1775 Double_t vertexZMax,
6acdbcb2 1776 Double_t ptTriggerMin,
1777 Double_t ptTriggerMax,
1778 Double_t ptAssociatedMin,
47e040b5 1779 Double_t ptAssociatedMax) {
a38fd7f3 1780 //Returns the 2D correlation function for (+-) pairs
622473b9 1781
1782 // security checks
b0d53e32 1783 if(psiMin > psiMax-0.00001){
1784 AliError("psiMax <= psiMin");
1785 return NULL;
1786 }
20006629 1787 if(vertexZMin > vertexZMax-0.00001){
1788 AliError("vertexZMax <= vertexZMin");
1789 return NULL;
1790 }
b0d53e32 1791 if(ptTriggerMin > ptTriggerMax-0.00001){
1792 AliError("ptTriggerMax <= ptTriggerMin");
1793 return NULL;
1794 }
1795 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1796 AliError("ptAssociatedMax <= ptAssociatedMin");
1797 return NULL;
1798 }
622473b9 1799
621329de 1800 // Psi_2: axis 0
622473b9 1801 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1802 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1803
20006629 1804 // Vz
1805 if(fVertexBinning){
1806 //Printf("Cutting on Vz...");
1807 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1808 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1809 }
1810
6acdbcb2 1811 // pt trigger
1812 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1813 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1814 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1815 }
621329de 1816
6acdbcb2 1817 // pt associated
1818 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1819 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1820
1821 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1822 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1823
1824 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1825 //TCanvas *c1 = new TCanvas("c1","");
1826 //c1->cd();
1827 //if(!gHistTest){
1828 //AliError("Projection of fHistP = NULL");
1829 //return gHistTest;
1830 //}
1831 //else{
1832 //gHistTest->DrawCopy("colz");
1833 //}
1834
1835 //0:step, 1: Delta eta, 2: Delta phi
621329de 1836 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1837 if(!gHist){
1838 AliError("Projection of fHistPN = NULL");
1839 return gHist;
1840 }
1841
6acdbcb2 1842 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1843 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1844 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1845
6acdbcb2 1846 //TCanvas *c2 = new TCanvas("c2","");
1847 //c2->cd();
1848 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1849
47e040b5 1850 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1851 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1852
1853 //normalize to bin width
1854 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1855
1856 return gHist;
1857}
1858
1859//____________________________________________________________________//
621329de 1860TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1861 Double_t psiMax,
20006629 1862 Double_t vertexZMin,
1863 Double_t vertexZMax,
6acdbcb2 1864 Double_t ptTriggerMin,
1865 Double_t ptTriggerMax,
1866 Double_t ptAssociatedMin,
47e040b5 1867 Double_t ptAssociatedMax) {
a38fd7f3 1868 //Returns the 2D correlation function for (+-) pairs
622473b9 1869
1870 // security checks
b0d53e32 1871 if(psiMin > psiMax-0.00001){
1872 AliError("psiMax <= psiMin");
1873 return NULL;
1874 }
20006629 1875 if(vertexZMin > vertexZMax-0.00001){
1876 AliError("vertexZMax <= vertexZMin");
1877 return NULL;
1878 }
b0d53e32 1879 if(ptTriggerMin > ptTriggerMax-0.00001){
1880 AliError("ptTriggerMax <= ptTriggerMin");
1881 return NULL;
1882 }
1883 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1884 AliError("ptAssociatedMax <= ptAssociatedMin");
1885 return NULL;
1886 }
622473b9 1887
621329de 1888 // Psi_2: axis 0
622473b9 1889 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1890 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 1891
1892 // Vz
1893 if(fVertexBinning){
1894 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1895 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1896 }
a38fd7f3 1897
6acdbcb2 1898 // pt trigger
1899 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1900 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1901 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1902 }
1903
1904 // pt associated
1905 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1906 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1907
1908 //0:step, 1: Delta eta, 2: Delta phi
621329de 1909 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1910 if(!gHist){
1911 AliError("Projection of fHistPN = NULL");
1912 return gHist;
1913 }
1914
a38fd7f3 1915 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1916 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 1917 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1918 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1919
1920 //normalize to bin width
1921 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1922
1923 return gHist;
1924}
1925
1926//____________________________________________________________________//
621329de 1927TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1928 Double_t psiMax,
20006629 1929 Double_t vertexZMin,
1930 Double_t vertexZMax,
6acdbcb2 1931 Double_t ptTriggerMin,
1932 Double_t ptTriggerMax,
1933 Double_t ptAssociatedMin,
47e040b5 1934 Double_t ptAssociatedMax) {
a38fd7f3 1935 //Returns the 2D correlation function for (+-) pairs
622473b9 1936
1937 // security checks
b0d53e32 1938 if(psiMin > psiMax-0.00001){
1939 AliError("psiMax <= psiMin");
1940 return NULL;
1941 }
20006629 1942 if(vertexZMin > vertexZMax-0.00001){
1943 AliError("vertexZMax <= vertexZMin");
1944 return NULL;
1945 }
b0d53e32 1946 if(ptTriggerMin > ptTriggerMax-0.00001){
1947 AliError("ptTriggerMax <= ptTriggerMin");
1948 return NULL;
1949 }
1950 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1951 AliError("ptAssociatedMax <= ptAssociatedMin");
1952 return NULL;
1953 }
622473b9 1954
621329de 1955 // Psi_2: axis 0
622473b9 1956 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1957 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1958
20006629 1959 // Vz
1960 if(fVertexBinning){
1961 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1962 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1963 }
1964
6acdbcb2 1965 // pt trigger
1966 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1967 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1968 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1969 }
1970
1971 // pt associated
1972 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1973 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1974
6acdbcb2 1975 //0:step, 1: Delta eta, 2: Delta phi
621329de 1976 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1977 if(!gHist){
1978 AliError("Projection of fHistPN = NULL");
1979 return gHist;
1980 }
1981
a38fd7f3 1982 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1983 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 1984 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1985 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1986
1987 //normalize to bin width
1988 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1989
1990 return gHist;
1991}
1992
1993//____________________________________________________________________//
621329de 1994TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1995 Double_t psiMax,
20006629 1996 Double_t vertexZMin,
1997 Double_t vertexZMax,
6acdbcb2 1998 Double_t ptTriggerMin,
1999 Double_t ptTriggerMax,
2000 Double_t ptAssociatedMin,
47e040b5 2001 Double_t ptAssociatedMax) {
a38fd7f3 2002 //Returns the 2D correlation function for (+-) pairs
622473b9 2003
2004 // security checks
2005 if(psiMin > psiMax-0.00001){
2006 AliError("psiMax <= psiMin");
2007 return NULL;
2008 }
20006629 2009 if(vertexZMin > vertexZMax-0.00001){
2010 AliError("vertexZMax <= vertexZMin");
2011 return NULL;
2012 }
622473b9 2013 if(ptTriggerMin > ptTriggerMax-0.00001){
2014 AliError("ptTriggerMax <= ptTriggerMin");
2015 return NULL;
2016 }
2017 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2018 AliError("ptAssociatedMax <= ptAssociatedMin");
2019 return NULL;
2020 }
2021
621329de 2022 // Psi_2: axis 0
622473b9 2023 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2024 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 2025
20006629 2026 // Vz
2027 if(fVertexBinning){
2028 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2029 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2030 }
2031
6acdbcb2 2032 // pt trigger
2033 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2034 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2035 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 2036 }
2037
2038 // pt associated
2039 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2040 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 2041
6acdbcb2 2042 //0:step, 1: Delta eta, 2: Delta phi
621329de 2043 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 2044 if(!gHist){
2045 AliError("Projection of fHistPN = NULL");
2046 return gHist;
2047 }
2048
a38fd7f3 2049 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2050 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 2051 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 2052 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 2053
2054 //normalize to bin width
2055 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 2056
2057 return gHist;
2058}
621329de 2059
7b610f27 2060//____________________________________________________________________//
2061TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2062 Double_t psiMax,
20006629 2063 Double_t vertexZMin,
2064 Double_t vertexZMax,
7b610f27 2065 Double_t ptTriggerMin,
2066 Double_t ptTriggerMax,
2067 Double_t ptAssociatedMin,
47e040b5 2068 Double_t ptAssociatedMax) {
7b610f27 2069 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 2070
2071 // security checks
b0d53e32 2072 if(psiMin > psiMax-0.00001){
2073 AliError("psiMax <= psiMin");
2074 return NULL;
2075 }
20006629 2076 if(vertexZMin > vertexZMax-0.00001){
2077 AliError("vertexZMax <= vertexZMin");
2078 return NULL;
2079 }
b0d53e32 2080 if(ptTriggerMin > ptTriggerMax-0.00001){
2081 AliError("ptTriggerMax <= ptTriggerMin");
2082 return NULL;
2083 }
2084 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2085 AliError("ptAssociatedMax <= ptAssociatedMin");
2086 return NULL;
2087 }
622473b9 2088
7b610f27 2089 // Psi_2: axis 0
622473b9 2090 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2091 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2092 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2093 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2094 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2095 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 2096
2097 // Vz
2098 if(fVertexBinning){
2099 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2100 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2101 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2102 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2103 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2104 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2105 }
7b610f27 2106
2107 // pt trigger
2108 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 2109 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2110 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2111 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2112 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2113 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2114 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 2115 }
2116
2117 // pt associated
2118 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 2119 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2120 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2121 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2122 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 2123
2124 //0:step, 1: Delta eta, 2: Delta phi
2125 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2126 if(!gHistNN){
2127 AliError("Projection of fHistNN = NULL");
2128 return gHistNN;
2129 }
2130 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2131 if(!gHistPP){
2132 AliError("Projection of fHistPP = NULL");
2133 return gHistPP;
2134 }
2135 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2136 if(!gHistNP){
2137 AliError("Projection of fHistNP = NULL");
2138 return gHistNP;
2139 }
2140 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2141 if(!gHistPN){
2142 AliError("Projection of fHistPN = NULL");
2143 return gHistPN;
2144 }
2145
2146 // sum all 2 particle histograms
2147 gHistNN->Add(gHistPP);
2148 gHistNN->Add(gHistNP);
2149 gHistNN->Add(gHistPN);
2150
47e040b5 2151 // divide by sum of + and - triggers
2152 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 2153 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 2154
2155 //normalize to bin width
2156 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 2157
2158 return gHistNN;
2159}
2160
f2e8af26 2161//____________________________________________________________________//
2162
4e023902 2163Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist,
f2e8af26 2164 Double_t &mean, Double_t &meanError,
2165 Double_t &sigma, Double_t &sigmaError,
2166 Double_t &skewness, Double_t &skewnessError,
2167 Double_t &kurtosis, Double_t &kurtosisError) {
2168 //
2169 // helper method to calculate the moments and errors of a TH1D anlytically
4e023902 2170 // fVariable = 1 for Delta eta (moments in full range)
2171 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
f2e8af26 2172 //
2173
2174 Bool_t success = kFALSE;
2175 mean = -1.;
2176 meanError = -1.;
2177 sigma = -1.;
2178 sigmaError = -1.;
2179 skewness = -1.;
2180 skewnessError = -1.;
2181 kurtosis = -1.;
2182 kurtosisError = -1.;
2183
2184 if(gHist){
2185
2186 // ----------------------------------------------------------------------
2187 // basic parameters of histogram
2188
2189 Int_t fNumberOfBins = gHist->GetNbinsX();
2190 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
2191
2192
2193 // ----------------------------------------------------------------------
2194 // first calculate the mean
2195
2196 Double_t fWeightedAverage = 0.;
f0754617 2197 Double_t fNormalization = 0.;
f2e8af26 2198
2199 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2200
4e023902 2201 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2202 if(fVariable == 2 &&
2203 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2204 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2205 continue;
2206 }
2207
f2e8af26 2208 fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
2209 fNormalization += gHist->GetBinContent(i);
f2e8af26 2210 }
2211
2212 mean = fWeightedAverage / fNormalization;
2213
f2e8af26 2214 // ----------------------------------------------------------------------
2215 // then calculate the higher moments
2216
a123fe62 2217 Double_t fMu = 0.;
2218 Double_t fMu2 = 0.;
2219 Double_t fMu3 = 0.;
2220 Double_t fMu4 = 0.;
2221 Double_t fMu5 = 0.;
2222 Double_t fMu6 = 0.;
2223 Double_t fMu7 = 0.;
2224 Double_t fMu8 = 0.;
f2e8af26 2225
2226 for(Int_t i = 1; i <= fNumberOfBins; i++) {
4e023902 2227
2228 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2229 if(fVariable == 2 &&
2230 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2231 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2232 continue;
2233 }
2234
a123fe62 2235 fMu += gHist->GetBinContent(i) * (gHist->GetBinCenter(i) - mean);
2236 fMu2 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2237 fMu3 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2238 fMu4 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2239 fMu5 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2240 fMu6 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2241 fMu7 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2242 fMu8 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
f2e8af26 2243 }
e47eeabd 2244
2245 // normalize to bin entries!
2246 fMu /= fNormalization;
2247 fMu2 /= fNormalization;
2248 fMu3 /= fNormalization;
2249 fMu4 /= fNormalization;
2250 fMu5 /= fNormalization;
2251 fMu6 /= fNormalization;
2252 fMu7 /= fNormalization;
2253 fMu8 /= fNormalization;
2254
2255 sigma = TMath::Sqrt(fMu2);
2256 skewness = fMu3 / TMath::Power(sigma,3);
2257 kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
f2e8af26 2258
f0754617 2259 // ----------------------------------------------------------------------
a123fe62 2260 // then calculate the higher moment errors
4e023902 2261 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2262 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2263 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
a123fe62 2264
2265 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2266
2267 // // assuming normal distribution (as done in ROOT)
2268 // meanError = sigma / TMath::Sqrt(normError);
2269 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2270 // skewnessError = TMath::Sqrt(6./(normError));
2271 // kurtosisError = TMath::Sqrt(24./(normError));
2272
2273 // use delta theorem paper (Luo - arXiv:1109.0593v1)
a9f20288 2274 Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2275 Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2276 Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
1af1cefb 2277 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2278 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2279 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
a123fe62 2280
4e023902 2281 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2282 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
a123fe62 2283
a9f20288 2284 if (TMath::Sqrt(normError) != 0){
2285 meanError = sigma / TMath::Sqrt(normError);
2286 }
2287 else return -999;
a123fe62 2288 sigmaError = TMath::Sqrt(Lambda11);
2289 skewnessError = TMath::Sqrt(Lambda22);
2290 kurtosisError = TMath::Sqrt(Lambda33);
a9f20288 2291
f0754617 2292
f2e8af26 2293 success = kTRUE;
2294 }
2295
2296
2297 return success;
2298}
2299
2300
73609a48 2301//____________________________________________________________________//
2302Float_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) {
2303 //
2304 // calculates dphistar
2305 //
2306 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2307
2308 static const Double_t kPi = TMath::Pi();
2309
2310 // circularity
2311// if (dphistar > 2 * kPi)
2312// dphistar -= 2 * kPi;
2313// if (dphistar < -2 * kPi)
2314// dphistar += 2 * kPi;
2315
2316 if (dphistar > kPi)
2317 dphistar = kPi * 2 - dphistar;
2318 if (dphistar < -kPi)
2319 dphistar = -kPi * 2 - dphistar;
2320 if (dphistar > kPi) // might look funny but is needed
2321 dphistar = kPi * 2 - dphistar;
2322
2323 return dphistar;
2324}
2325
2922bbe3 2326
2327