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