]> 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
updated
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalancePsi.cxx
CommitLineData
0879e280 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14/* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
15
16//-----------------------------------------------------------------
17// Balance Function class
18// This is the class to deal with the Balance Function wrt Psi analysis
19// Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21
22
23//ROOT
24#include <Riostream.h>
621329de 25#include <TCanvas.h>
0879e280 26#include <TMath.h>
27#include <TAxis.h>
28#include <TH2D.h>
29#include <TH3D.h>
30#include <TLorentzVector.h>
31#include <TObjArray.h>
32#include <TGraphErrors.h>
33#include <TString.h>
34
35#include "AliVParticle.h"
36#include "AliMCParticle.h"
37#include "AliESDtrack.h"
38#include "AliAODTrack.h"
9afe3098 39#include "AliTHn.h"
35aff0f3 40#include "AliAnalysisTaskTriggeredBF.h"
0879e280 41
42#include "AliBalancePsi.h"
43
44ClassImp(AliBalancePsi)
45
46//____________________________________________________________________//
47AliBalancePsi::AliBalancePsi() :
48 TObject(),
9fd4b54e 49 fShuffle(kFALSE),
0879e280 50 fAnalysisLevel("ESD"),
51 fAnalyzedEvents(0) ,
52 fCentralityId(0) ,
53 fCentStart(0.),
54 fCentStop(0.),
9afe3098 55 fHistP(0),
56 fHistN(0),
57 fHistPN(0),
58 fHistNP(0),
59 fHistPP(0),
60 fHistNN(0),
73609a48 61 fHistHBTbefore(0),
62 fHistHBTafter(0),
63 fHistConversionbefore(0),
64 fHistConversionafter(0),
c8ad9635 65 fHistPsiMinusPhi(0),
f0754617 66 fHistResonancesBefore(0),
67 fHistResonancesRho(0),
68 fHistResonancesK0(0),
69 fHistResonancesLambda(0),
5a861dc6 70 fHistQbefore(0),
71 fHistQafter(0),
73609a48 72 fPsiInterval(15.),
7c04a4d2 73 fDeltaEtaMax(2.0),
5a861dc6 74 fResonancesCut(kFALSE),
75 fHBTCut(kFALSE),
76 fConversionCut(kFALSE),
77 fQCut(kFALSE),
683cbfb5 78 fVertexBinning(kFALSE),
f0fb4ac1 79 fEventClass("EventPlane"){
0879e280 80 // Default constructor
0879e280 81}
82
0879e280 83//____________________________________________________________________//
84AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
9fd4b54e 85 TObject(balance), fShuffle(balance.fShuffle),
0879e280 86 fAnalysisLevel(balance.fAnalysisLevel),
87 fAnalyzedEvents(balance.fAnalyzedEvents),
88 fCentralityId(balance.fCentralityId),
89 fCentStart(balance.fCentStart),
90 fCentStop(balance.fCentStop),
9afe3098 91 fHistP(balance.fHistP),
92 fHistN(balance.fHistN),
93 fHistPN(balance.fHistPN),
94 fHistNP(balance.fHistNP),
95 fHistPP(balance.fHistPP),
96 fHistNN(balance.fHistNN),
73609a48 97 fHistHBTbefore(balance.fHistHBTbefore),
98 fHistHBTafter(balance.fHistHBTafter),
99 fHistConversionbefore(balance.fHistConversionbefore),
100 fHistConversionafter(balance.fHistConversionafter),
c8ad9635 101 fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
f0754617 102 fHistResonancesBefore(balance.fHistResonancesBefore),
103 fHistResonancesRho(balance.fHistResonancesRho),
104 fHistResonancesK0(balance.fHistResonancesK0),
105 fHistResonancesLambda(balance.fHistResonancesLambda),
5a861dc6 106 fHistQbefore(balance.fHistQbefore),
107 fHistQafter(balance.fHistQafter),
73609a48 108 fPsiInterval(balance.fPsiInterval),
7c04a4d2 109 fDeltaEtaMax(balance.fDeltaEtaMax),
2922bbe3 110 fResonancesCut(balance.fResonancesCut),
73609a48 111 fHBTCut(balance.fHBTCut),
f0fb4ac1 112 fConversionCut(balance.fConversionCut),
5a861dc6 113 fQCut(balance.fQCut),
683cbfb5 114 fVertexBinning(balance.fVertexBinning),
f0fb4ac1 115 fEventClass("EventPlane"){
0879e280 116 //copy constructor
0879e280 117}
118
119//____________________________________________________________________//
120AliBalancePsi::~AliBalancePsi() {
121 // Destructor
9afe3098 122 delete fHistP;
123 delete fHistN;
124 delete fHistPN;
125 delete fHistNP;
126 delete fHistPP;
127 delete fHistNN;
73609a48 128
129 delete fHistHBTbefore;
130 delete fHistHBTafter;
131 delete fHistConversionbefore;
132 delete fHistConversionafter;
c8ad9635 133 delete fHistPsiMinusPhi;
f0754617 134 delete fHistResonancesBefore;
135 delete fHistResonancesRho;
136 delete fHistResonancesK0;
137 delete fHistResonancesLambda;
5a861dc6 138 delete fHistQbefore;
139 delete fHistQafter;
f0fb4ac1 140
0879e280 141}
142
143//____________________________________________________________________//
9afe3098 144void AliBalancePsi::InitHistograms() {
145 // single particle histograms
211b716d 146
147 // global switch disabling the reference
148 // (to avoid "Replacing existing TH1" if several wagons are created in train)
149 Bool_t oldStatus = TH1::AddDirectoryStatus();
150 TH1::AddDirectory(kFALSE);
151
9afe3098 152 Int_t anaSteps = 1; // analysis steps
9fd4b54e 153 Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
154 Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
155 TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
9afe3098 156
157 // two particle histograms
9fd4b54e 158 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
159 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
160 TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
f0fb4ac1 161 /**********************************************************
162
163 ======> Modification: Change Event Classification Scheme
164
165 ---> fEventClass == "EventPlane"
166
167 Default operation with Event Plane
168
169 ---> fEventClass == "Multiplicity"
170
171 Work with kTPCITStracklet multiplicity (from GetReferenceMultiplicity)
172
173 ---> fEventClass == "Centrality"
174
175 Work with Centrality Bins
176
177 ***********************************************************/
178
179 //--- Multiplicity Bins ------------------------------------
180 const Int_t kMultBins = 8;
181 //A first rough attempt at four bins
182 Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80};
183 //----------------------------------------------------------
184
185 //--- Centrality Bins --------------------------------------
683cbfb5 186 const Int_t kNCentralityBins = 9;
187 const Int_t kNCentralityBinsVertex = 26;
188 Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
189 Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.};
f0fb4ac1 190 //----------------------------------------------------------
191
192 //--- Event Plane Bins -------------------------------------
193 //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest)
194 const Int_t kNPsi2Bins = 4;
195 Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
196 //----------------------------------------------------------
197
198 //Depending on fEventClass Variable, do one thing or the other...
199 if(fEventClass == "Multiplicity"){
200 iBinSingle[0] = kMultBins;
201 dBinsSingle[0] = kMultBinLimits;
202 axisTitleSingle[0] = "kTPCITStracklet multiplicity";
203 iBinPair[0] = kMultBins;
204 dBinsPair[0] = kMultBinLimits;
205 axisTitlePair[0] = "kTPCITStracklet multiplicity";
206 }
207 if(fEventClass == "Centrality"){
683cbfb5 208 // fine binning in case of vertex Z binning
209 if(fVertexBinning){
210 iBinSingle[0] = kNCentralityBinsVertex;
211 dBinsSingle[0] = centralityBinsVertex;
212
213 iBinPair[0] = kNCentralityBinsVertex;
214 dBinsPair[0] = centralityBinsVertex;
215 }
216 else{
217 iBinSingle[0] = kNCentralityBins;
218 dBinsSingle[0] = centralityBins;
219
220 iBinPair[0] = kNCentralityBins;
221 dBinsPair[0] = centralityBins;
222 }
f0fb4ac1 223 axisTitleSingle[0] = "Centrality percentile [%]";
f0fb4ac1 224 axisTitlePair[0] = "Centrality percentile [%]";
225 }
226 if(fEventClass == "EventPlane"){
227 iBinSingle[0] = kNPsi2Bins;
228 dBinsSingle[0] = psi2Bins;
229 axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
230 iBinPair[0] = kNPsi2Bins;
231 dBinsPair[0] = psi2Bins;
232 axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
233 }
9afe3098 234
683cbfb5 235 // delta eta
236 const Int_t kNDeltaEtaBins = 80;
237 const Int_t kNDeltaEtaBinsVertex = 40;
238 Double_t deltaEtaBins[kNDeltaEtaBins+1];
239 Double_t deltaEtaBinsVertex[kNDeltaEtaBinsVertex+1];
240 for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
241 deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins;
242 for(Int_t i = 0; i < kNDeltaEtaBinsVertex+1; i++)
243 deltaEtaBinsVertex[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBinsVertex;
244
245 // coarse binning in case of vertex Z binning
246 if(fVertexBinning){
247 iBinPair[1] = kNDeltaEtaBinsVertex;
248 dBinsPair[1] = deltaEtaBinsVertex;
249 }
250 else{
251 iBinPair[1] = kNDeltaEtaBins;
252 dBinsPair[1] = deltaEtaBins;
253 }
254 axisTitlePair[1] = "#Delta#eta";
621329de 255
256 // delta phi
9afe3098 257 const Int_t kNDeltaPhiBins = 72;
258 Double_t deltaPhiBins[kNDeltaPhiBins+1];
259 for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
c8ad9635 260 //deltaPhiBins[i] = -180.0 + i * 5.;
261 deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
9afe3098 262 }
621329de 263 iBinPair[2] = kNDeltaPhiBins;
264 dBinsPair[2] = deltaPhiBins;
c8ad9635 265 axisTitlePair[2] = "#Delta#varphi (rad)";
9afe3098 266
a38fd7f3 267 // pt(trigger-associated)
683cbfb5 268 const Int_t kNPtBins = 16;
269 const Int_t kNPtBinsVertex = 5;
270 Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
271 Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0};
272
273 // coarse binning in case of vertex Z binning
274 if(fVertexBinning){
275 iBinSingle[1] = kNPtBinsVertex;
276 dBinsSingle[1] = ptBinsVertex;
277
278 iBinPair[3] = kNPtBinsVertex;
279 dBinsPair[3] = ptBinsVertex;
280
281 iBinPair[4] = kNPtBinsVertex;
282 dBinsPair[4] = ptBinsVertex;
283 }
284 else{
285 iBinSingle[1] = kNPtBins;
286 dBinsSingle[1] = ptBins;
287
288 iBinPair[3] = kNPtBins;
289 dBinsPair[3] = ptBins;
290
291 iBinPair[4] = kNPtBins;
292 dBinsPair[4] = ptBins;
293 }
294
c8ad9635 295 axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
683cbfb5 296 axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
297 axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
298
299 // vertex Z
300 const Int_t kNVertexZBins = 1;
301 const Int_t kNVertexZBinsVertex = 9;
302 Double_t vertexZBins[kNVertexZBins+1] = {-10., 10.};
303 Double_t vertexZBinsVertex[kNVertexZBinsVertex+1] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.};
304
305 // vertex Z binning or not
306 if(fVertexBinning){
307 iBinSingle[2] = kNVertexZBinsVertex;
308 dBinsSingle[2] = vertexZBinsVertex;
309
310 iBinPair[5] = kNVertexZBinsVertex;
311 dBinsPair[5] = vertexZBinsVertex;
312 }
313 else{
314 iBinSingle[2] = kNVertexZBins;
315 dBinsSingle[2] = vertexZBins;
621329de 316
683cbfb5 317 iBinPair[5] = kNVertexZBins;
318 dBinsPair[5] = vertexZBins;
319 }
320
321 axisTitleSingle[2] = "v_{Z} (cm)";
322 axisTitlePair[5] = "v_{Z} (cm)";
a38fd7f3 323
9afe3098 324
325 TString histName;
326 //+ triggered particles
327 histName = "fHistP";
9fd4b54e 328 if(fShuffle) histName.Append("_shuffle");
9afe3098 329 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 330 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
331 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 332 fHistP->SetBinLimits(j, dBinsSingle[j]);
333 fHistP->SetVarTitle(j, axisTitleSingle[j]);
0879e280 334 }
9afe3098 335
336 //- triggered particles
337 histName = "fHistN";
9fd4b54e 338 if(fShuffle) histName.Append("_shuffle");
9afe3098 339 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 340 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
341 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 342 fHistN->SetBinLimits(j, dBinsSingle[j]);
343 fHistN->SetVarTitle(j, axisTitleSingle[j]);
0879e280 344 }
9afe3098 345
346 //+- pairs
347 histName = "fHistPN";
9fd4b54e 348 if(fShuffle) histName.Append("_shuffle");
9afe3098 349 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 350 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
351 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 352 fHistPN->SetBinLimits(j, dBinsPair[j]);
353 fHistPN->SetVarTitle(j, axisTitlePair[j]);
0879e280 354 }
0879e280 355
9afe3098 356 //-+ pairs
357 histName = "fHistNP";
9fd4b54e 358 if(fShuffle) histName.Append("_shuffle");
9afe3098 359 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 360 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
361 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 362 fHistNP->SetBinLimits(j, dBinsPair[j]);
363 fHistNP->SetVarTitle(j, axisTitlePair[j]);
0879e280 364 }
0879e280 365
9afe3098 366 //++ pairs
367 histName = "fHistPP";
9fd4b54e 368 if(fShuffle) histName.Append("_shuffle");
9afe3098 369 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 370 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
371 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 372 fHistPP->SetBinLimits(j, dBinsPair[j]);
373 fHistPP->SetVarTitle(j, axisTitlePair[j]);
374 }
375
376 //-- pairs
377 histName = "fHistNN";
9fd4b54e 378 if(fShuffle) histName.Append("_shuffle");
9afe3098 379 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 380 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
381 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 382 fHistNN->SetBinLimits(j, dBinsPair[j]);
383 fHistNN->SetVarTitle(j, axisTitlePair[j]);
0879e280 384 }
f06d59b3 385 AliInfo("Finished setting up the AliTHn");
73609a48 386
387 // QA histograms
c8ad9635 388 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
389 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
390 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
391 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
f0754617 392 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
393 fHistResonancesBefore = new TH3D("fHistResonancesBefore","before resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
394 fHistResonancesRho = new TH3D("fHistResonancesRho","after #rho resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
395 fHistResonancesK0 = new TH3D("fHistResonancesK0","after #rho, K0 resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
396 fHistResonancesLambda = new TH3D("fHistResonancesLambda","after #rho, K0, Lambda resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
5a861dc6 397 fHistQbefore = new TH3D("fHistQbefore","before momentum difference cut;#Delta#eta;#Delta#phi;q = p_{1} - p_{2} (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;q = p_{1} - p_{2} (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
668 Double_t qMin = 0.1; //const for the time being (should be changeable later on)
669 TVector3 vDiff,v[2];
670 Double_t momentumDifference = 0.;
671
672 v[0].SetPtEtaPhi(firstPt,firstEta,firstPhi);
673 v[1].SetPtEtaPhi(secondPt[j],secondEta[j],secondPhi[j]);
674 vDiff = v[1] - v[0];
675 momentumDifference = TMath::Abs(vDiff.Mag());
676
677 fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],momentumDifference);
678 if(momentumDifference < qMin) continue;
679 fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],momentumDifference);
680
681 }
682
35aff0f3 683 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
684 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
685 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
686 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
8795e993 687 else {
688 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
689 continue;
690 }
6acdbcb2 691 }//end of 2nd particle loop
692 }//end of 1st particle loop
0879e280 693}
694
0879e280 695//____________________________________________________________________//
9afe3098 696TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
697 Int_t iVariablePair,
9afe3098 698 Double_t psiMin,
6acdbcb2 699 Double_t psiMax,
20006629 700 Double_t vertexZMin,
701 Double_t vertexZMax,
6acdbcb2 702 Double_t ptTriggerMin,
703 Double_t ptTriggerMax,
704 Double_t ptAssociatedMin,
705 Double_t ptAssociatedMax) {
9afe3098 706 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 707 //(private members) of the AliBalancePsi class.
621329de 708 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
709 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 710
711 // security checks
b0d53e32 712 if(psiMin > psiMax-0.00001){
713 AliError("psiMax <= psiMin");
714 return NULL;
715 }
20006629 716 if(vertexZMin > vertexZMax-0.00001){
717 AliError("vertexZMax <= vertexZMin");
718 return NULL;
719 }
b0d53e32 720 if(ptTriggerMin > ptTriggerMax-0.00001){
721 AliError("ptTriggerMax <= ptTriggerMin");
722 return NULL;
723 }
724 if(ptAssociatedMin > ptAssociatedMax-0.00001){
725 AliError("ptAssociatedMax <= ptAssociatedMin");
726 return NULL;
727 }
622473b9 728
9afe3098 729 // Psi_2
622473b9 730 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
731 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
732 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
733 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
734 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
735 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
9afe3098 736
20006629 737 // Vz
738 if(fVertexBinning){
739 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
740 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
741 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
742 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
743 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
744 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
745 }
746
6acdbcb2 747 // pt trigger
748 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 749 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
750 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
751 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
752 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
753 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
754 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 755 }
756
757 // pt associated
758 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 759 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
760 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
761 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
762 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 763 }
764
a38fd7f3 765 //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));
766
9afe3098 767 // Project into the wanted space (1st: analysis step, 2nd: axis)
20006629 768 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); //
769 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); //
770 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); //
771 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); //
772 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); //
773 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); //
9afe3098 774
775 TH1D *gHistBalanceFunctionHistogram = 0x0;
776 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
777 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
778 gHistBalanceFunctionHistogram->Reset();
779
780 switch(iVariablePair) {
621329de 781 case 1:
c8ad9635 782 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
783 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 784 break;
621329de 785 case 2:
c8ad9635 786 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
787 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 788 break;
9afe3098 789 default:
790 break;
791 }
0879e280 792
0879e280 793 hTemp1->Sumw2();
794 hTemp2->Sumw2();
795 hTemp3->Sumw2();
796 hTemp4->Sumw2();
a38fd7f3 797 hTemp1->Add(hTemp3,-1.);
0879e280 798 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 799 hTemp2->Add(hTemp4,-1.);
0879e280 800 hTemp2->Scale(1./hTemp6->GetEntries());
801 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 802 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 803
804 //normalize to bin width
805 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
0879e280 806 }
807
0879e280 808 return gHistBalanceFunctionHistogram;
809}
a38fd7f3 810
f0c5040c 811//____________________________________________________________________//
812TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
813 Int_t iVariablePair,
814 Double_t psiMin,
815 Double_t psiMax,
20006629 816 Double_t vertexZMin,
817 Double_t vertexZMax,
f0c5040c 818 Double_t ptTriggerMin,
819 Double_t ptTriggerMax,
820 Double_t ptAssociatedMin,
821 Double_t ptAssociatedMax,
822 AliBalancePsi *bfMix) {
823 //Returns the BF histogram, extracted from the 6 AliTHn objects
824 //after dividing each correlation function by the Event Mixing one
825 //(private members) of the AliBalancePsi class.
826 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
827 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 828
829 // security checks
b0d53e32 830 if(psiMin > psiMax-0.00001){
831 AliError("psiMax <= psiMin");
832 return NULL;
833 }
20006629 834 if(vertexZMin > vertexZMax-0.00001){
835 AliError("vertexZMax <= vertexZMin");
836 return NULL;
837 }
b0d53e32 838 if(ptTriggerMin > ptTriggerMax-0.00001){
839 AliError("ptTriggerMax <= ptTriggerMin");
840 return NULL;
841 }
842 if(ptAssociatedMin > ptAssociatedMax-0.00001){
843 AliError("ptAssociatedMax <= ptAssociatedMin");
844 return NULL;
845 }
622473b9 846
f0c5040c 847 // Psi_2
622473b9 848 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
849 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
850 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
851 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
852 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
853 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 854
20006629 855 // Vz
856 if(fVertexBinning){
857 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
858 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
859 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
860 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
861 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
862 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
863 }
864
f0c5040c 865 // pt trigger
866 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 867 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
868 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
869 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
870 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
871 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
872 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 873 }
874
875 // pt associated
876 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 877 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
878 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
879 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
880 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 881 }
882
883 // ============================================================================================
884 // the same for event mixing
885 AliTHn *fHistPMix = bfMix->GetHistNp();
886 AliTHn *fHistNMix = bfMix->GetHistNn();
887 AliTHn *fHistPNMix = bfMix->GetHistNpn();
888 AliTHn *fHistNPMix = bfMix->GetHistNnp();
889 AliTHn *fHistPPMix = bfMix->GetHistNpp();
890 AliTHn *fHistNNMix = bfMix->GetHistNnn();
891
892 // Psi_2
622473b9 893 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
894 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
895 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
896 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
897 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
898 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 899
20006629 900 // Vz
901 if(fVertexBinning){
902 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
903 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
904 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
905 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
906 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
907 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
908 }
909
f0c5040c 910 // pt trigger
911 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 912 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
913 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
914 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
915 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
916 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
917 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 918 }
919
920 // pt associated
921 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 922 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
923 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
924 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
925 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 926 }
927 // ============================================================================================
928
929 //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));
930
931 // Project into the wanted space (1st: analysis step, 2nd: axis)
932 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
933 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
934 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
935 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
936 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
937 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
938
939 // ============================================================================================
940 // the same for event mixing
f0fb4ac1 941 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
942 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
943 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
944 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
945 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
946 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
f0c5040c 947 // ============================================================================================
948
949 TH1D *gHistBalanceFunctionHistogram = 0x0;
950 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
951 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
952 gHistBalanceFunctionHistogram->Reset();
953
954 switch(iVariablePair) {
955 case 1:
956 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
957 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
958 break;
959 case 2:
960 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
961 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
962 break;
963 default:
964 break;
965 }
966
967 hTemp1->Sumw2();
968 hTemp2->Sumw2();
969 hTemp3->Sumw2();
970 hTemp4->Sumw2();
971 hTemp1Mix->Sumw2();
972 hTemp2Mix->Sumw2();
973 hTemp3Mix->Sumw2();
974 hTemp4Mix->Sumw2();
975
976 hTemp1->Scale(1./hTemp5->GetEntries());
977 hTemp3->Scale(1./hTemp5->GetEntries());
978 hTemp2->Scale(1./hTemp6->GetEntries());
979 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 980 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
981 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
982 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
983 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 984
985 hTemp1->Divide(hTemp1Mix);
986 hTemp2->Divide(hTemp2Mix);
987 hTemp3->Divide(hTemp3Mix);
988 hTemp4->Divide(hTemp4Mix);
989
990 hTemp1->Add(hTemp3,-1.);
991 hTemp2->Add(hTemp4,-1.);
992
993 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
994 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 995
996 //normalize to bin width
997 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 998 }
999
1000 return gHistBalanceFunctionHistogram;
1001}
1002
a38fd7f3 1003//____________________________________________________________________//
621329de 1004TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 1005 Double_t psiMax,
20006629 1006 Double_t vertexZMin,
1007 Double_t vertexZMax,
6acdbcb2 1008 Double_t ptTriggerMin,
1009 Double_t ptTriggerMax,
1010 Double_t ptAssociatedMin,
1011 Double_t ptAssociatedMax) {
621329de 1012 //Returns the BF histogram in Delta eta vs Delta phi,
1013 //extracted from the 6 AliTHn objects
1014 //(private members) of the AliBalancePsi class.
1015 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1016 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 1017
1018 // security checks
b0d53e32 1019 if(psiMin > psiMax-0.00001){
1020 AliError("psiMax <= psiMin");
1021 return NULL;
1022 }
20006629 1023 if(vertexZMin > vertexZMax-0.00001){
1024 AliError("vertexZMax <= vertexZMin");
1025 return NULL;
1026 }
b0d53e32 1027 if(ptTriggerMin > ptTriggerMax-0.00001){
1028 AliError("ptTriggerMax <= ptTriggerMin");
1029 return NULL;
1030 }
1031 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1032 AliError("ptAssociatedMax <= ptAssociatedMin");
1033 return NULL;
1034 }
622473b9 1035
621329de 1036 TString histName = "gHistBalanceFunctionHistogram2D";
1037
1038 // Psi_2
622473b9 1039 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1040 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1041 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1042 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1043 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1044 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1045
20006629 1046 // Vz
1047 if(fVertexBinning){
1048 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1049 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1050 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1051 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1052 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1053 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1054 }
1055
6acdbcb2 1056 // pt trigger
1057 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1058 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1059 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1060 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1061 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1062 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1063 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1064 }
1065
1066 // pt associated
1067 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1068 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1069 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1070 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1071 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1072 }
1073
1074 //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 1075
1076 // Project into the wanted space (1st: analysis step, 2nd: axis)
1077 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1078 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1079 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1080 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1081 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1082 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1083
1084 TH2D *gHistBalanceFunctionHistogram = 0x0;
1085 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1086 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1087 gHistBalanceFunctionHistogram->Reset();
c8ad9635 1088 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1089 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1090 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 1091
1092 hTemp1->Sumw2();
1093 hTemp2->Sumw2();
1094 hTemp3->Sumw2();
1095 hTemp4->Sumw2();
1096 hTemp1->Add(hTemp3,-1.);
1097 hTemp1->Scale(1./hTemp5->GetEntries());
1098 hTemp2->Add(hTemp4,-1.);
1099 hTemp2->Scale(1./hTemp6->GetEntries());
1100 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1101 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1102
1103 //normalize to bin width
1104 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
621329de 1105 }
1106
1107 return gHistBalanceFunctionHistogram;
1108}
1109
f0c5040c 1110//____________________________________________________________________//
1111TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1112 Double_t psiMax,
20006629 1113 Double_t vertexZMin,
1114 Double_t vertexZMax,
f0c5040c 1115 Double_t ptTriggerMin,
1116 Double_t ptTriggerMax,
1117 Double_t ptAssociatedMin,
1118 Double_t ptAssociatedMax,
1119 AliBalancePsi *bfMix) {
1120 //Returns the BF histogram in Delta eta vs Delta phi,
1121 //after dividing each correlation function by the Event Mixing one
1122 //extracted from the 6 AliTHn objects
1123 //(private members) of the AliBalancePsi class.
1124 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1125 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1126 TString histName = "gHistBalanceFunctionHistogram2D";
1127
1128 if(!bfMix){
1129 AliError("balance function object for event mixing not available");
1130 return NULL;
1131 }
1132
622473b9 1133 // security checks
b0d53e32 1134 if(psiMin > psiMax-0.00001){
1135 AliError("psiMax <= psiMin");
1136 return NULL;
1137 }
20006629 1138 if(vertexZMin > vertexZMax-0.00001){
1139 AliError("vertexZMax <= vertexZMin");
1140 return NULL;
1141 }
b0d53e32 1142 if(ptTriggerMin > ptTriggerMax-0.00001){
1143 AliError("ptTriggerMax <= ptTriggerMin");
1144 return NULL;
1145 }
1146 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1147 AliError("ptAssociatedMax <= ptAssociatedMin");
1148 return NULL;
1149 }
622473b9 1150
f0c5040c 1151 // Psi_2
622473b9 1152 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1153 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1154 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1155 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1156 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1157 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 1158
1159 // Vz
1160 if(fVertexBinning){
1161 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1162 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1163 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1164 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1165 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1166 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1167 }
f0c5040c 1168
1169 // pt trigger
1170 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1171 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1172 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1173 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1174 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1175 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1176 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1177 }
1178
1179 // pt associated
1180 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1181 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1182 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1183 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1184 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1185 }
1186
1187
1188 // ============================================================================================
1189 // the same for event mixing
1190 AliTHn *fHistPMix = bfMix->GetHistNp();
1191 AliTHn *fHistNMix = bfMix->GetHistNn();
1192 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1193 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1194 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1195 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1196
1197 // Psi_2
622473b9 1198 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1199 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1200 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1201 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1202 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1203 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1204
20006629 1205 // Vz
1206 if(fVertexBinning){
1207 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1208 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1209 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1210 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1211 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1212 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1213 }
f0c5040c 1214 // pt trigger
1215 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1216 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1217 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1218 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1219 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1220 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1221 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1222 }
1223
1224 // pt associated
1225 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1226 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1227 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1228 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1229 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1230 }
1231 // ============================================================================================
1232
1233
1234 //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)));
1235
1236 // Project into the wanted space (1st: analysis step, 2nd: axis)
1237 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1238 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1239 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1240 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1241 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1242 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1243
1244 // ============================================================================================
1245 // the same for event mixing
1246 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1247 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1248 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1249 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
f0fb4ac1 1250 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1251 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
f0c5040c 1252 // ============================================================================================
1253
1254 TH2D *gHistBalanceFunctionHistogram = 0x0;
1255 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1256 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1257 gHistBalanceFunctionHistogram->Reset();
1258 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1259 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1260 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1261
1262 hTemp1->Sumw2();
1263 hTemp2->Sumw2();
1264 hTemp3->Sumw2();
1265 hTemp4->Sumw2();
1266 hTemp1Mix->Sumw2();
1267 hTemp2Mix->Sumw2();
1268 hTemp3Mix->Sumw2();
1269 hTemp4Mix->Sumw2();
1270
1271 hTemp1->Scale(1./hTemp5->GetEntries());
1272 hTemp3->Scale(1./hTemp5->GetEntries());
1273 hTemp2->Scale(1./hTemp6->GetEntries());
1274 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 1275 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1276 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1277 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1278 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 1279
1280 hTemp1->Divide(hTemp1Mix);
1281 hTemp2->Divide(hTemp2Mix);
1282 hTemp3->Divide(hTemp3Mix);
1283 hTemp4->Divide(hTemp4Mix);
1284
1285 hTemp1->Add(hTemp3,-1.);
1286 hTemp2->Add(hTemp4,-1.);
1287
1288 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1289 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1290
1291 //normalize to bin width
1292 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 1293 }
1294
1295 return gHistBalanceFunctionHistogram;
1296}
1297
f2e8af26 1298//____________________________________________________________________//
1299TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1300 Double_t psiMin,
1301 Double_t psiMax,
20006629 1302 Double_t vertexZMin,
1303 Double_t vertexZMax,
f2e8af26 1304 Double_t ptTriggerMin,
1305 Double_t ptTriggerMax,
1306 Double_t ptAssociatedMin,
1307 Double_t ptAssociatedMax,
1308 AliBalancePsi *bfMix) {
1309 //Returns the BF histogram in Delta eta OR Delta phi,
1310 //after dividing each correlation function by the Event Mixing one
1311 // (But the division is done here in 2D, this was basically done to check the Integral)
1312 //extracted from the 6 AliTHn objects
1313 //(private members) of the AliBalancePsi class.
1314 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1315 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1316 TString histName = "gHistBalanceFunctionHistogram1D";
1317
1318 if(!bfMix){
1319 AliError("balance function object for event mixing not available");
1320 return NULL;
1321 }
1322
622473b9 1323 // security checks
b0d53e32 1324 if(psiMin > psiMax-0.00001){
1325 AliError("psiMax <= psiMin");
1326 return NULL;
1327 }
20006629 1328 if(vertexZMin > vertexZMax-0.00001){
1329 AliError("vertexZMax <= vertexZMin");
1330 return NULL;
1331 }
b0d53e32 1332 if(ptTriggerMin > ptTriggerMax-0.00001){
1333 AliError("ptTriggerMax <= ptTriggerMin");
1334 return NULL;
1335 }
1336 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1337 AliError("ptAssociatedMax <= ptAssociatedMin");
1338 return NULL;
1339 }
622473b9 1340
f2e8af26 1341 // Psi_2
622473b9 1342 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1343 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1344 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1345 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1346 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1347 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1348
20006629 1349 // Vz
1350 if(fVertexBinning){
1351 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1352 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1353 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1354 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1355 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1356 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1357 }
1358
f2e8af26 1359 // pt trigger
1360 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1361 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1362 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1363 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1364 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1365 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1366 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1367 }
1368
1369 // pt associated
1370 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1371 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1372 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1373 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1374 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1375 }
1376
1377
1378 // ============================================================================================
1379 // the same for event mixing
1380 AliTHn *fHistPMix = bfMix->GetHistNp();
1381 AliTHn *fHistNMix = bfMix->GetHistNn();
1382 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1383 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1384 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1385 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1386
1387 // Psi_2
622473b9 1388 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1389 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1390 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1391 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1392 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1393 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1394
20006629 1395 // Vz
1396 if(fVertexBinning){
1397 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1398 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1399 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1400 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1401 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1402 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1403 }
1404
f2e8af26 1405 // pt trigger
1406 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1407 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1408 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1409 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1410 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1411 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1412 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1413 }
1414
1415 // pt associated
1416 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1417 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1418 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1419 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1420 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1421 }
1422 // ============================================================================================
1423
1424
1425 //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)));
1426
1427 // Project into the wanted space (1st: analysis step, 2nd: axis)
1428 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1429 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1430 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1431 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1432 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1433 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1434
1435 // ============================================================================================
1436 // the same for event mixing
1437 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1438 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1439 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1440 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1441 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1442 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1443 // ============================================================================================
1444
1445 TH1D *gHistBalanceFunctionHistogram = 0x0;
1446 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1447
1448 hTemp1->Sumw2();
1449 hTemp2->Sumw2();
1450 hTemp3->Sumw2();
1451 hTemp4->Sumw2();
1452 hTemp1Mix->Sumw2();
1453 hTemp2Mix->Sumw2();
1454 hTemp3Mix->Sumw2();
1455 hTemp4Mix->Sumw2();
1456
1457 hTemp1->Scale(1./hTemp5->GetEntries());
1458 hTemp3->Scale(1./hTemp5->GetEntries());
1459 hTemp2->Scale(1./hTemp6->GetEntries());
1460 hTemp4->Scale(1./hTemp6->GetEntries());
1461
1462 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1463 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1464 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1465 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
1466
1467 hTemp1->Divide(hTemp1Mix);
1468 hTemp2->Divide(hTemp2Mix);
1469 hTemp3->Divide(hTemp3Mix);
1470 hTemp4->Divide(hTemp4Mix);
1471
1472 // now only project on one axis
1473 TH1D *h1DTemp1 = NULL;
1474 TH1D *h1DTemp2 = NULL;
1475 TH1D *h1DTemp3 = NULL;
1476 TH1D *h1DTemp4 = NULL;
1477 if(!bPhi){
1478 h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
1479 h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
1480 h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
1481 h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
1482 }
1483 else{
1484 h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
1485 h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
1486 h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
1487 h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
1488 }
1489
1490 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1491 gHistBalanceFunctionHistogram->Reset();
1492 if(!bPhi){
1493 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1494 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1495 }
1496 else{
1497 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1498 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1499 }
1500
1501 h1DTemp1->Add(h1DTemp3,-1.);
1502 h1DTemp2->Add(h1DTemp4,-1.);
1503
1504 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1505 gHistBalanceFunctionHistogram->Scale(0.5);
1506 }
1507
1508 return gHistBalanceFunctionHistogram;
1509}
1510
f0c5040c 1511
621329de 1512//____________________________________________________________________//
1513TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1514 Double_t psiMax,
20006629 1515 Double_t vertexZMin,
1516 Double_t vertexZMax,
6acdbcb2 1517 Double_t ptTriggerMin,
1518 Double_t ptTriggerMax,
1519 Double_t ptAssociatedMin,
47e040b5 1520 Double_t ptAssociatedMax) {
a38fd7f3 1521 //Returns the 2D correlation function for (+-) pairs
622473b9 1522
1523 // security checks
b0d53e32 1524 if(psiMin > psiMax-0.00001){
1525 AliError("psiMax <= psiMin");
1526 return NULL;
1527 }
20006629 1528 if(vertexZMin > vertexZMax-0.00001){
1529 AliError("vertexZMax <= vertexZMin");
1530 return NULL;
1531 }
b0d53e32 1532 if(ptTriggerMin > ptTriggerMax-0.00001){
1533 AliError("ptTriggerMax <= ptTriggerMin");
1534 return NULL;
1535 }
1536 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1537 AliError("ptAssociatedMax <= ptAssociatedMin");
1538 return NULL;
1539 }
622473b9 1540
621329de 1541 // Psi_2: axis 0
622473b9 1542 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1543 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1544
20006629 1545 // Vz
1546 if(fVertexBinning){
1547 //Printf("Cutting on Vz...");
1548 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1549 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1550 }
1551
6acdbcb2 1552 // pt trigger
1553 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1554 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1555 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1556 }
621329de 1557
6acdbcb2 1558 // pt associated
1559 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1560 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1561
1562 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1563 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1564
1565 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1566 //TCanvas *c1 = new TCanvas("c1","");
1567 //c1->cd();
1568 //if(!gHistTest){
1569 //AliError("Projection of fHistP = NULL");
1570 //return gHistTest;
1571 //}
1572 //else{
1573 //gHistTest->DrawCopy("colz");
1574 //}
1575
1576 //0:step, 1: Delta eta, 2: Delta phi
621329de 1577 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1578 if(!gHist){
1579 AliError("Projection of fHistPN = NULL");
1580 return gHist;
1581 }
1582
6acdbcb2 1583 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1584 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1585 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1586
6acdbcb2 1587 //TCanvas *c2 = new TCanvas("c2","");
1588 //c2->cd();
1589 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1590
47e040b5 1591 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1592 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1593
1594 //normalize to bin width
1595 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1596
1597 return gHist;
1598}
1599
1600//____________________________________________________________________//
621329de 1601TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1602 Double_t psiMax,
20006629 1603 Double_t vertexZMin,
1604 Double_t vertexZMax,
6acdbcb2 1605 Double_t ptTriggerMin,
1606 Double_t ptTriggerMax,
1607 Double_t ptAssociatedMin,
47e040b5 1608 Double_t ptAssociatedMax) {
a38fd7f3 1609 //Returns the 2D correlation function for (+-) pairs
622473b9 1610
1611 // security checks
b0d53e32 1612 if(psiMin > psiMax-0.00001){
1613 AliError("psiMax <= psiMin");
1614 return NULL;
1615 }
20006629 1616 if(vertexZMin > vertexZMax-0.00001){
1617 AliError("vertexZMax <= vertexZMin");
1618 return NULL;
1619 }
b0d53e32 1620 if(ptTriggerMin > ptTriggerMax-0.00001){
1621 AliError("ptTriggerMax <= ptTriggerMin");
1622 return NULL;
1623 }
1624 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1625 AliError("ptAssociatedMax <= ptAssociatedMin");
1626 return NULL;
1627 }
622473b9 1628
621329de 1629 // Psi_2: axis 0
622473b9 1630 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1631 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 1632
1633 // Vz
1634 if(fVertexBinning){
1635 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1636 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1637 }
a38fd7f3 1638
6acdbcb2 1639 // pt trigger
1640 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1641 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1642 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1643 }
1644
1645 // pt associated
1646 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1647 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1648
1649 //0:step, 1: Delta eta, 2: Delta phi
621329de 1650 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1651 if(!gHist){
1652 AliError("Projection of fHistPN = NULL");
1653 return gHist;
1654 }
1655
a38fd7f3 1656 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1657 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 1658 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1659 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1660
1661 //normalize to bin width
1662 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1663
1664 return gHist;
1665}
1666
1667//____________________________________________________________________//
621329de 1668TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1669 Double_t psiMax,
20006629 1670 Double_t vertexZMin,
1671 Double_t vertexZMax,
6acdbcb2 1672 Double_t ptTriggerMin,
1673 Double_t ptTriggerMax,
1674 Double_t ptAssociatedMin,
47e040b5 1675 Double_t ptAssociatedMax) {
a38fd7f3 1676 //Returns the 2D correlation function for (+-) pairs
622473b9 1677
1678 // security checks
b0d53e32 1679 if(psiMin > psiMax-0.00001){
1680 AliError("psiMax <= psiMin");
1681 return NULL;
1682 }
20006629 1683 if(vertexZMin > vertexZMax-0.00001){
1684 AliError("vertexZMax <= vertexZMin");
1685 return NULL;
1686 }
b0d53e32 1687 if(ptTriggerMin > ptTriggerMax-0.00001){
1688 AliError("ptTriggerMax <= ptTriggerMin");
1689 return NULL;
1690 }
1691 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1692 AliError("ptAssociatedMax <= ptAssociatedMin");
1693 return NULL;
1694 }
622473b9 1695
621329de 1696 // Psi_2: axis 0
622473b9 1697 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1698 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1699
20006629 1700 // Vz
1701 if(fVertexBinning){
1702 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1703 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1704 }
1705
6acdbcb2 1706 // pt trigger
1707 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1708 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1709 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1710 }
1711
1712 // pt associated
1713 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1714 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1715
6acdbcb2 1716 //0:step, 1: Delta eta, 2: Delta phi
621329de 1717 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1718 if(!gHist){
1719 AliError("Projection of fHistPN = NULL");
1720 return gHist;
1721 }
1722
a38fd7f3 1723 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1724 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 1725 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1726 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1727
1728 //normalize to bin width
1729 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1730
1731 return gHist;
1732}
1733
1734//____________________________________________________________________//
621329de 1735TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1736 Double_t psiMax,
20006629 1737 Double_t vertexZMin,
1738 Double_t vertexZMax,
6acdbcb2 1739 Double_t ptTriggerMin,
1740 Double_t ptTriggerMax,
1741 Double_t ptAssociatedMin,
47e040b5 1742 Double_t ptAssociatedMax) {
a38fd7f3 1743 //Returns the 2D correlation function for (+-) pairs
622473b9 1744
1745 // security checks
1746 if(psiMin > psiMax-0.00001){
1747 AliError("psiMax <= psiMin");
1748 return NULL;
1749 }
20006629 1750 if(vertexZMin > vertexZMax-0.00001){
1751 AliError("vertexZMax <= vertexZMin");
1752 return NULL;
1753 }
622473b9 1754 if(ptTriggerMin > ptTriggerMax-0.00001){
1755 AliError("ptTriggerMax <= ptTriggerMin");
1756 return NULL;
1757 }
1758 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1759 AliError("ptAssociatedMax <= ptAssociatedMin");
1760 return NULL;
1761 }
1762
621329de 1763 // Psi_2: axis 0
622473b9 1764 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1765 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1766
20006629 1767 // Vz
1768 if(fVertexBinning){
1769 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1770 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1771 }
1772
6acdbcb2 1773 // pt trigger
1774 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1775 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1776 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1777 }
1778
1779 // pt associated
1780 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1781 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1782
6acdbcb2 1783 //0:step, 1: Delta eta, 2: Delta phi
621329de 1784 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 1785 if(!gHist){
1786 AliError("Projection of fHistPN = NULL");
1787 return gHist;
1788 }
1789
a38fd7f3 1790 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1791 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 1792 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1793 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1794
1795 //normalize to bin width
1796 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1797
1798 return gHist;
1799}
621329de 1800
7b610f27 1801//____________________________________________________________________//
1802TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1803 Double_t psiMax,
20006629 1804 Double_t vertexZMin,
1805 Double_t vertexZMax,
7b610f27 1806 Double_t ptTriggerMin,
1807 Double_t ptTriggerMax,
1808 Double_t ptAssociatedMin,
47e040b5 1809 Double_t ptAssociatedMax) {
7b610f27 1810 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 1811
1812 // security checks
b0d53e32 1813 if(psiMin > psiMax-0.00001){
1814 AliError("psiMax <= psiMin");
1815 return NULL;
1816 }
20006629 1817 if(vertexZMin > vertexZMax-0.00001){
1818 AliError("vertexZMax <= vertexZMin");
1819 return NULL;
1820 }
b0d53e32 1821 if(ptTriggerMin > ptTriggerMax-0.00001){
1822 AliError("ptTriggerMax <= ptTriggerMin");
1823 return NULL;
1824 }
1825 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1826 AliError("ptAssociatedMax <= ptAssociatedMin");
1827 return NULL;
1828 }
622473b9 1829
7b610f27 1830 // Psi_2: axis 0
622473b9 1831 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1832 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1833 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1834 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1835 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1836 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
20006629 1837
1838 // Vz
1839 if(fVertexBinning){
1840 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1841 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1842 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1843 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1844 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1845 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1846 }
7b610f27 1847
1848 // pt trigger
1849 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1850 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1851 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1852 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1853 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1854 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1855 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 1856 }
1857
1858 // pt associated
1859 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1860 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1861 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1862 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1863 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 1864
1865 //0:step, 1: Delta eta, 2: Delta phi
1866 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1867 if(!gHistNN){
1868 AliError("Projection of fHistNN = NULL");
1869 return gHistNN;
1870 }
1871 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1872 if(!gHistPP){
1873 AliError("Projection of fHistPP = NULL");
1874 return gHistPP;
1875 }
1876 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1877 if(!gHistNP){
1878 AliError("Projection of fHistNP = NULL");
1879 return gHistNP;
1880 }
1881 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1882 if(!gHistPN){
1883 AliError("Projection of fHistPN = NULL");
1884 return gHistPN;
1885 }
1886
1887 // sum all 2 particle histograms
1888 gHistNN->Add(gHistPP);
1889 gHistNN->Add(gHistNP);
1890 gHistNN->Add(gHistPN);
1891
47e040b5 1892 // divide by sum of + and - triggers
1893 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 1894 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 1895
1896 //normalize to bin width
1897 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 1898
1899 return gHistNN;
1900}
1901
f2e8af26 1902//____________________________________________________________________//
1903
4e023902 1904Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist,
f2e8af26 1905 Double_t &mean, Double_t &meanError,
1906 Double_t &sigma, Double_t &sigmaError,
1907 Double_t &skewness, Double_t &skewnessError,
1908 Double_t &kurtosis, Double_t &kurtosisError) {
1909 //
1910 // helper method to calculate the moments and errors of a TH1D anlytically
4e023902 1911 // fVariable = 1 for Delta eta (moments in full range)
1912 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
f2e8af26 1913 //
1914
1915 Bool_t success = kFALSE;
1916 mean = -1.;
1917 meanError = -1.;
1918 sigma = -1.;
1919 sigmaError = -1.;
1920 skewness = -1.;
1921 skewnessError = -1.;
1922 kurtosis = -1.;
1923 kurtosisError = -1.;
1924
1925 if(gHist){
1926
1927 // ----------------------------------------------------------------------
1928 // basic parameters of histogram
1929
1930 Int_t fNumberOfBins = gHist->GetNbinsX();
1931 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
1932
1933
1934 // ----------------------------------------------------------------------
1935 // first calculate the mean
1936
1937 Double_t fWeightedAverage = 0.;
f0754617 1938 Double_t fNormalization = 0.;
f2e8af26 1939
1940 for(Int_t i = 1; i <= fNumberOfBins; i++) {
1941
4e023902 1942 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
1943 if(fVariable == 2 &&
1944 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
1945 gHist->GetBinCenter(i) > TMath::Pi()/2)){
1946 continue;
1947 }
1948
f2e8af26 1949 fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
1950 fNormalization += gHist->GetBinContent(i);
f2e8af26 1951 }
1952
1953 mean = fWeightedAverage / fNormalization;
1954
f2e8af26 1955 // ----------------------------------------------------------------------
1956 // then calculate the higher moments
1957
a123fe62 1958 Double_t fMu = 0.;
1959 Double_t fMu2 = 0.;
1960 Double_t fMu3 = 0.;
1961 Double_t fMu4 = 0.;
1962 Double_t fMu5 = 0.;
1963 Double_t fMu6 = 0.;
1964 Double_t fMu7 = 0.;
1965 Double_t fMu8 = 0.;
f2e8af26 1966
1967 for(Int_t i = 1; i <= fNumberOfBins; i++) {
4e023902 1968
1969 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
1970 if(fVariable == 2 &&
1971 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
1972 gHist->GetBinCenter(i) > TMath::Pi()/2)){
1973 continue;
1974 }
1975
a123fe62 1976 fMu += gHist->GetBinContent(i) * (gHist->GetBinCenter(i) - mean);
1977 fMu2 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
1978 fMu3 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
1979 fMu4 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
1980 fMu5 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
1981 fMu6 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
1982 fMu7 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
1983 fMu8 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
f2e8af26 1984 }
1985
a123fe62 1986 sigma = TMath::Sqrt(fMu2 / fNormalization);
1987 skewness = fMu3 / fNormalization / TMath::Power(sigma,3);
1988 kurtosis = fMu4 / fNormalization / TMath::Power(sigma,4) - 3;
f2e8af26 1989
f0754617 1990 // ----------------------------------------------------------------------
a123fe62 1991 // then calculate the higher moment errors
4e023902 1992 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
1993 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
1994 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
a123fe62 1995
1996 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
1997
1998 // // assuming normal distribution (as done in ROOT)
1999 // meanError = sigma / TMath::Sqrt(normError);
2000 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2001 // skewnessError = TMath::Sqrt(6./(normError));
2002 // kurtosisError = TMath::Sqrt(24./(normError));
2003
2004 // use delta theorem paper (Luo - arXiv:1109.0593v1)
2005 Double_t Lambda11 = (fMu4-1)*sigma*sigma/(4*normError);
2006 Double_t Lambda22 = (9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError;
2007 Double_t Lambda33 = (-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError;
1af1cefb 2008 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2009 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2010 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
a123fe62 2011
4e023902 2012 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2013 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
a123fe62 2014
2015 meanError = sigma / TMath::Sqrt(normError);
2016 sigmaError = TMath::Sqrt(Lambda11);
2017 skewnessError = TMath::Sqrt(Lambda22);
2018 kurtosisError = TMath::Sqrt(Lambda33);
2019
f0754617 2020
f2e8af26 2021 success = kTRUE;
2022 }
2023
2024
2025 return success;
2026}
2027
2028
73609a48 2029//____________________________________________________________________//
2030Float_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) {
2031 //
2032 // calculates dphistar
2033 //
2034 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2035
2036 static const Double_t kPi = TMath::Pi();
2037
2038 // circularity
2039// if (dphistar > 2 * kPi)
2040// dphistar -= 2 * kPi;
2041// if (dphistar < -2 * kPi)
2042// dphistar += 2 * kPi;
2043
2044 if (dphistar > kPi)
2045 dphistar = kPi * 2 - dphistar;
2046 if (dphistar < -kPi)
2047 dphistar = -kPi * 2 - dphistar;
2048 if (dphistar > kPi) // might look funny but is needed
2049 dphistar = kPi * 2 - dphistar;
2050
2051 return dphistar;
2052}
2053
2922bbe3 2054
2055