]> 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
QA histograms for resonance cuts
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalancePsi.cxx
CommitLineData
0879e280 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14/* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
15
16//-----------------------------------------------------------------
17// Balance Function class
18// This is the class to deal with the Balance Function wrt Psi analysis
19// Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21
22
23//ROOT
24#include <Riostream.h>
621329de 25#include <TCanvas.h>
0879e280 26#include <TMath.h>
27#include <TAxis.h>
28#include <TH2D.h>
29#include <TH3D.h>
30#include <TLorentzVector.h>
31#include <TObjArray.h>
32#include <TGraphErrors.h>
33#include <TString.h>
34
35#include "AliVParticle.h"
36#include "AliMCParticle.h"
37#include "AliESDtrack.h"
38#include "AliAODTrack.h"
9afe3098 39#include "AliTHn.h"
35aff0f3 40#include "AliAnalysisTaskTriggeredBF.h"
0879e280 41
42#include "AliBalancePsi.h"
43
44ClassImp(AliBalancePsi)
45
46//____________________________________________________________________//
47AliBalancePsi::AliBalancePsi() :
48 TObject(),
9fd4b54e 49 fShuffle(kFALSE),
0879e280 50 fAnalysisLevel("ESD"),
51 fAnalyzedEvents(0) ,
52 fCentralityId(0) ,
53 fCentStart(0.),
54 fCentStop(0.),
9afe3098 55 fHistP(0),
56 fHistN(0),
57 fHistPN(0),
58 fHistNP(0),
59 fHistPP(0),
60 fHistNN(0),
73609a48 61 fHistHBTbefore(0),
62 fHistHBTafter(0),
63 fHistConversionbefore(0),
64 fHistConversionafter(0),
c8ad9635 65 fHistPsiMinusPhi(0),
f0754617 66 fHistResonancesBefore(0),
67 fHistResonancesRho(0),
68 fHistResonancesK0(0),
69 fHistResonancesLambda(0),
73609a48 70 fPsiInterval(15.),
7c04a4d2 71 fDeltaEtaMax(2.0),
2922bbe3 72 fResonancesCut(kTRUE),
73 fHBTCut(kTRUE),
74 fConversionCut(kTRUE),
683cbfb5 75 fVertexBinning(kFALSE),
f0fb4ac1 76 fEventClass("EventPlane"){
0879e280 77 // Default constructor
0879e280 78}
79
0879e280 80//____________________________________________________________________//
81AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
9fd4b54e 82 TObject(balance), fShuffle(balance.fShuffle),
0879e280 83 fAnalysisLevel(balance.fAnalysisLevel),
84 fAnalyzedEvents(balance.fAnalyzedEvents),
85 fCentralityId(balance.fCentralityId),
86 fCentStart(balance.fCentStart),
87 fCentStop(balance.fCentStop),
9afe3098 88 fHistP(balance.fHistP),
89 fHistN(balance.fHistN),
90 fHistPN(balance.fHistPN),
91 fHistNP(balance.fHistNP),
92 fHistPP(balance.fHistPP),
93 fHistNN(balance.fHistNN),
73609a48 94 fHistHBTbefore(balance.fHistHBTbefore),
95 fHistHBTafter(balance.fHistHBTafter),
96 fHistConversionbefore(balance.fHistConversionbefore),
97 fHistConversionafter(balance.fHistConversionafter),
c8ad9635 98 fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
f0754617 99 fHistResonancesBefore(balance.fHistResonancesBefore),
100 fHistResonancesRho(balance.fHistResonancesRho),
101 fHistResonancesK0(balance.fHistResonancesK0),
102 fHistResonancesLambda(balance.fHistResonancesLambda),
73609a48 103 fPsiInterval(balance.fPsiInterval),
7c04a4d2 104 fDeltaEtaMax(balance.fDeltaEtaMax),
2922bbe3 105 fResonancesCut(balance.fResonancesCut),
73609a48 106 fHBTCut(balance.fHBTCut),
f0fb4ac1 107 fConversionCut(balance.fConversionCut),
683cbfb5 108 fVertexBinning(balance.fVertexBinning),
f0fb4ac1 109 fEventClass("EventPlane"){
0879e280 110 //copy constructor
0879e280 111}
112
113//____________________________________________________________________//
114AliBalancePsi::~AliBalancePsi() {
115 // Destructor
9afe3098 116 delete fHistP;
117 delete fHistN;
118 delete fHistPN;
119 delete fHistNP;
120 delete fHistPP;
121 delete fHistNN;
73609a48 122
123 delete fHistHBTbefore;
124 delete fHistHBTafter;
125 delete fHistConversionbefore;
126 delete fHistConversionafter;
c8ad9635 127 delete fHistPsiMinusPhi;
f0754617 128 delete fHistResonancesBefore;
129 delete fHistResonancesRho;
130 delete fHistResonancesK0;
131 delete fHistResonancesLambda;
f0fb4ac1 132
0879e280 133}
134
135//____________________________________________________________________//
9afe3098 136void AliBalancePsi::InitHistograms() {
137 // single particle histograms
211b716d 138
139 // global switch disabling the reference
140 // (to avoid "Replacing existing TH1" if several wagons are created in train)
141 Bool_t oldStatus = TH1::AddDirectoryStatus();
142 TH1::AddDirectory(kFALSE);
143
9afe3098 144 Int_t anaSteps = 1; // analysis steps
9fd4b54e 145 Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
146 Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
147 TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
9afe3098 148
149 // two particle histograms
9fd4b54e 150 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
151 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
152 TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
f0fb4ac1 153 /**********************************************************
154
155 ======> Modification: Change Event Classification Scheme
156
157 ---> fEventClass == "EventPlane"
158
159 Default operation with Event Plane
160
161 ---> fEventClass == "Multiplicity"
162
163 Work with kTPCITStracklet multiplicity (from GetReferenceMultiplicity)
164
165 ---> fEventClass == "Centrality"
166
167 Work with Centrality Bins
168
169 ***********************************************************/
170
171 //--- Multiplicity Bins ------------------------------------
172 const Int_t kMultBins = 8;
173 //A first rough attempt at four bins
174 Double_t kMultBinLimits[kMultBins+1]={0,10,20,30,40,50,60,70,80};
175 //----------------------------------------------------------
176
177 //--- Centrality Bins --------------------------------------
683cbfb5 178 const Int_t kNCentralityBins = 9;
179 const Int_t kNCentralityBinsVertex = 26;
180 Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
181 Double_t centralityBinsVertex[kNCentralityBinsVertex+1] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.};
f0fb4ac1 182 //----------------------------------------------------------
183
184 //--- Event Plane Bins -------------------------------------
185 //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)
186 const Int_t kNPsi2Bins = 4;
187 Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
188 //----------------------------------------------------------
189
190 //Depending on fEventClass Variable, do one thing or the other...
191 if(fEventClass == "Multiplicity"){
192 iBinSingle[0] = kMultBins;
193 dBinsSingle[0] = kMultBinLimits;
194 axisTitleSingle[0] = "kTPCITStracklet multiplicity";
195 iBinPair[0] = kMultBins;
196 dBinsPair[0] = kMultBinLimits;
197 axisTitlePair[0] = "kTPCITStracklet multiplicity";
198 }
199 if(fEventClass == "Centrality"){
683cbfb5 200 // fine binning in case of vertex Z binning
201 if(fVertexBinning){
202 iBinSingle[0] = kNCentralityBinsVertex;
203 dBinsSingle[0] = centralityBinsVertex;
204
205 iBinPair[0] = kNCentralityBinsVertex;
206 dBinsPair[0] = centralityBinsVertex;
207 }
208 else{
209 iBinSingle[0] = kNCentralityBins;
210 dBinsSingle[0] = centralityBins;
211
212 iBinPair[0] = kNCentralityBins;
213 dBinsPair[0] = centralityBins;
214 }
f0fb4ac1 215 axisTitleSingle[0] = "Centrality percentile [%]";
f0fb4ac1 216 axisTitlePair[0] = "Centrality percentile [%]";
217 }
218 if(fEventClass == "EventPlane"){
219 iBinSingle[0] = kNPsi2Bins;
220 dBinsSingle[0] = psi2Bins;
221 axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
222 iBinPair[0] = kNPsi2Bins;
223 dBinsPair[0] = psi2Bins;
224 axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
225 }
9afe3098 226
683cbfb5 227 // delta eta
228 const Int_t kNDeltaEtaBins = 80;
229 const Int_t kNDeltaEtaBinsVertex = 40;
230 Double_t deltaEtaBins[kNDeltaEtaBins+1];
231 Double_t deltaEtaBinsVertex[kNDeltaEtaBinsVertex+1];
232 for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
233 deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins;
234 for(Int_t i = 0; i < kNDeltaEtaBinsVertex+1; i++)
235 deltaEtaBinsVertex[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBinsVertex;
236
237 // coarse binning in case of vertex Z binning
238 if(fVertexBinning){
239 iBinPair[1] = kNDeltaEtaBinsVertex;
240 dBinsPair[1] = deltaEtaBinsVertex;
241 }
242 else{
243 iBinPair[1] = kNDeltaEtaBins;
244 dBinsPair[1] = deltaEtaBins;
245 }
246 axisTitlePair[1] = "#Delta#eta";
621329de 247
248 // delta phi
9afe3098 249 const Int_t kNDeltaPhiBins = 72;
250 Double_t deltaPhiBins[kNDeltaPhiBins+1];
251 for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
c8ad9635 252 //deltaPhiBins[i] = -180.0 + i * 5.;
253 deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
9afe3098 254 }
621329de 255 iBinPair[2] = kNDeltaPhiBins;
256 dBinsPair[2] = deltaPhiBins;
c8ad9635 257 axisTitlePair[2] = "#Delta#varphi (rad)";
9afe3098 258
a38fd7f3 259 // pt(trigger-associated)
683cbfb5 260 const Int_t kNPtBins = 16;
261 const Int_t kNPtBinsVertex = 5;
262 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.};
263 Double_t ptBinsVertex[kNPtBinsVertex+1] = {0.2,1.0,2.0,3.0,4.0,8.0};
264
265 // coarse binning in case of vertex Z binning
266 if(fVertexBinning){
267 iBinSingle[1] = kNPtBinsVertex;
268 dBinsSingle[1] = ptBinsVertex;
269
270 iBinPair[3] = kNPtBinsVertex;
271 dBinsPair[3] = ptBinsVertex;
272
273 iBinPair[4] = kNPtBinsVertex;
274 dBinsPair[4] = ptBinsVertex;
275 }
276 else{
277 iBinSingle[1] = kNPtBins;
278 dBinsSingle[1] = ptBins;
279
280 iBinPair[3] = kNPtBins;
281 dBinsPair[3] = ptBins;
282
283 iBinPair[4] = kNPtBins;
284 dBinsPair[4] = ptBins;
285 }
286
c8ad9635 287 axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
683cbfb5 288 axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
289 axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
290
291 // vertex Z
292 const Int_t kNVertexZBins = 1;
293 const Int_t kNVertexZBinsVertex = 9;
294 Double_t vertexZBins[kNVertexZBins+1] = {-10., 10.};
295 Double_t vertexZBinsVertex[kNVertexZBinsVertex+1] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.};
296
297 // vertex Z binning or not
298 if(fVertexBinning){
299 iBinSingle[2] = kNVertexZBinsVertex;
300 dBinsSingle[2] = vertexZBinsVertex;
301
302 iBinPair[5] = kNVertexZBinsVertex;
303 dBinsPair[5] = vertexZBinsVertex;
304 }
305 else{
306 iBinSingle[2] = kNVertexZBins;
307 dBinsSingle[2] = vertexZBins;
621329de 308
683cbfb5 309 iBinPair[5] = kNVertexZBins;
310 dBinsPair[5] = vertexZBins;
311 }
312
313 axisTitleSingle[2] = "v_{Z} (cm)";
314 axisTitlePair[5] = "v_{Z} (cm)";
a38fd7f3 315
9afe3098 316
317 TString histName;
318 //+ triggered particles
319 histName = "fHistP";
9fd4b54e 320 if(fShuffle) histName.Append("_shuffle");
9afe3098 321 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 322 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
323 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 324 fHistP->SetBinLimits(j, dBinsSingle[j]);
325 fHistP->SetVarTitle(j, axisTitleSingle[j]);
0879e280 326 }
9afe3098 327
328 //- triggered particles
329 histName = "fHistN";
9fd4b54e 330 if(fShuffle) histName.Append("_shuffle");
9afe3098 331 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 332 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
333 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 334 fHistN->SetBinLimits(j, dBinsSingle[j]);
335 fHistN->SetVarTitle(j, axisTitleSingle[j]);
0879e280 336 }
9afe3098 337
338 //+- pairs
339 histName = "fHistPN";
9fd4b54e 340 if(fShuffle) histName.Append("_shuffle");
9afe3098 341 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 342 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
343 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 344 fHistPN->SetBinLimits(j, dBinsPair[j]);
345 fHistPN->SetVarTitle(j, axisTitlePair[j]);
0879e280 346 }
0879e280 347
9afe3098 348 //-+ pairs
349 histName = "fHistNP";
9fd4b54e 350 if(fShuffle) histName.Append("_shuffle");
9afe3098 351 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 352 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
353 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 354 fHistNP->SetBinLimits(j, dBinsPair[j]);
355 fHistNP->SetVarTitle(j, axisTitlePair[j]);
0879e280 356 }
0879e280 357
9afe3098 358 //++ pairs
359 histName = "fHistPP";
9fd4b54e 360 if(fShuffle) histName.Append("_shuffle");
9afe3098 361 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 362 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
363 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 364 fHistPP->SetBinLimits(j, dBinsPair[j]);
365 fHistPP->SetVarTitle(j, axisTitlePair[j]);
366 }
367
368 //-- pairs
369 histName = "fHistNN";
9fd4b54e 370 if(fShuffle) histName.Append("_shuffle");
9afe3098 371 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 372 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
373 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 374 fHistNN->SetBinLimits(j, dBinsPair[j]);
375 fHistNN->SetVarTitle(j, axisTitlePair[j]);
0879e280 376 }
f06d59b3 377 AliInfo("Finished setting up the AliTHn");
73609a48 378
379 // QA histograms
c8ad9635 380 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
381 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
382 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
383 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
f0754617 384 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
385 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);
386 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);
387 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);
388 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);
389
211b716d 390
391 TH1::AddDirectory(oldStatus);
392
0879e280 393}
394
395//____________________________________________________________________//
f06d59b3 396void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
397 TObjArray *particles,
73609a48 398 TObjArray *particlesMixed,
683cbfb5 399 Float_t bSign,
400 Double_t kMultorCent,
401 Double_t vertexZ) {
0879e280 402 // Calculates the balance function
403 fAnalyzedEvents++;
f0fb4ac1 404
0879e280 405 // Initialize histograms if not done yet
9afe3098 406 if(!fHistPN){
0879e280 407 AliWarning("Histograms not yet initialized! --> Will be done now");
408 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
409 InitHistograms();
410 }
411
9fd4b54e 412 Double_t trackVariablesSingle[kTrackVariablesSingle];
413 Double_t trackVariablesPair[kTrackVariablesPair];
f06d59b3 414
415 if (!particles){
416 AliWarning("particles TObjArray is NULL pointer --> return");
417 return;
418 }
9afe3098 419
f06d59b3 420 // define end of particle loops
421 Int_t iMax = particles->GetEntriesFast();
422 Int_t jMax = iMax;
423 if (particlesMixed)
424 jMax = particlesMixed->GetEntriesFast();
425
426 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
427 TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
428
429 TArrayF secondEta(jMax);
430 TArrayF secondPhi(jMax);
431 TArrayF secondPt(jMax);
c1847400 432 TArrayS secondCharge(jMax);
35aff0f3 433 TArrayD secondCorrection(jMax);
f06d59b3 434
435 for (Int_t i=0; i<jMax; i++){
436 secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
437 secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
438 secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
c1847400 439 secondCharge[i] = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
35aff0f3 440 secondCorrection[i] = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction(); //==========================correction
f06d59b3 441 }
442
2922bbe3 443 //TLorenzVector implementation for resonances
444 TLorentzVector vectorMother, vectorDaughter[2];
445 TParticle pPion, pProton, pRho0, pK0s, pLambda;
446 pPion.SetPdgCode(211); //pion
447 pRho0.SetPdgCode(113); //rho0
448 pK0s.SetPdgCode(310); //K0s
449 pProton.SetPdgCode(2212); //proton
450 pLambda.SetPdgCode(3122); //Lambda
451 Double_t gWidthForRho0 = 0.01;
452 Double_t gWidthForK0s = 0.01;
453 Double_t gWidthForLambda = 0.006;
454 Double_t nSigmaRejection = 3.0;
455
f06d59b3 456 // 1st particle loop
73609a48 457 for (Int_t i = 0; i < iMax; i++) {
35aff0f3 458 //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
459 AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
6acdbcb2 460
461 // some optimization
462 Float_t firstEta = firstParticle->Eta();
463 Float_t firstPhi = firstParticle->Phi();
464 Float_t firstPt = firstParticle->Pt();
35aff0f3 465 Float_t firstCorrection = firstParticle->Correction();//==========================correction
466
6acdbcb2 467 // Event plane (determine psi bin)
468 Double_t gPsiMinusPhi = 0.;
469 Double_t gPsiMinusPhiBin = -10.;
470 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
471 //in-plane
c8ad9635 472 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
473 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
6acdbcb2 474 gPsiMinusPhiBin = 0.0;
475 //intermediate
c8ad9635 476 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
477 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
478 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
479 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
6acdbcb2 480 gPsiMinusPhiBin = 1.0;
481 //out of plane
c8ad9635 482 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
483 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
6acdbcb2 484 gPsiMinusPhiBin = 2.0;
485 //everything else
486 else
487 gPsiMinusPhiBin = 3.0;
488
c8ad9635 489 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
73609a48 490
491 Short_t charge1 = (Short_t) firstParticle->Charge();
6acdbcb2 492
493 trackVariablesSingle[0] = gPsiMinusPhiBin;
f0fb4ac1 494 trackVariablesSingle[1] = firstPt;
bb473a33 495 if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
683cbfb5 496 trackVariablesSingle[2] = vertexZ;
497
6acdbcb2 498
499 //fill single particle histograms
35aff0f3 500 if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
501 else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
6acdbcb2 502
c1847400 503 // 2nd particle loop
504 for(Int_t j = 0; j < jMax; j++) {
505
506 if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
507
d26c6047 508 // pT,Assoc < pT,Trig
509 if(firstPt < secondPt[j]) continue;
510
c1847400 511 Short_t charge2 = secondCharge[j];
9afe3098 512
f0fb4ac1 513 trackVariablesPair[0] = trackVariablesSingle[0];
6acdbcb2 514 trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
515 trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
c8ad9635 516 //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
517 //trackVariablesPair[2] -= 360.;
518 //if (trackVariablesPair[2] < - 180.)
519 //trackVariablesPair[2] += 360.;
520 if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
521 trackVariablesPair[2] -= 2.*TMath::Pi();
522 if (trackVariablesPair[2] < - TMath::Pi())
523 trackVariablesPair[2] += 2.*TMath::Pi();
524 if (trackVariablesPair[2] < - TMath::Pi()/2.)
525 trackVariablesPair[2] += 2.*TMath::Pi();
9afe3098 526
6acdbcb2 527 trackVariablesPair[3] = firstPt; // pt trigger
528 trackVariablesPair[4] = secondPt[j]; // pt
683cbfb5 529 trackVariablesPair[5] = vertexZ; // z of the primary vertex
2922bbe3 530
531 //Exclude resonances for the calculation of pairs by looking
532 //at the invariant mass and not considering the pairs that
533 //fall within 3sigma from the mass peak of: rho0, K0s, Lambda
534 if(fResonancesCut) {
f0754617 535
2922bbe3 536 //rho0
537 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
538 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
539 vectorMother = vectorDaughter[0] + vectorDaughter[1];
f0754617 540 fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
2922bbe3 541 if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0)
542 continue;
f0754617 543 fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
2922bbe3 544
545 //K0s
546 if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s)
547 continue;
f0754617 548 fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
549
2922bbe3 550
551 //Lambda
552 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
553 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass());
554 vectorMother = vectorDaughter[0] + vectorDaughter[1];
555 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
556 continue;
557
558 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass());
559 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
560 vectorMother = vectorDaughter[0] + vectorDaughter[1];
561 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
562 continue;
f0754617 563 fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
564
2922bbe3 565 }//resonance cut
73609a48 566
567 // HBT like cut
da8b2d30 568 if(fHBTCut){ // VERSION 3 (all pairs)
569 //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
73609a48 570 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
571 // continue;
572
573 Double_t deta = firstEta - secondEta[j];
574 Double_t dphi = firstPhi - secondPhi[j];
575 // VERSION 2 (Taken from DPhiCorrelations)
576 // the variables & cuthave been developed by the HBT group
577 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
578 fHistHBTbefore->Fill(deta,dphi);
579
580 // optimization
581 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
582 {
583 // phi in rad
c8ad9635 584 //Float_t phi1rad = firstPhi*TMath::DegToRad();
585 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
586 Float_t phi1rad = firstPhi;
587 Float_t phi2rad = secondPhi[j];
73609a48 588
589 // check first boundaries to see if is worth to loop and find the minimum
590 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
591 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
592
593 const Float_t kLimit = 0.02 * 3;
594
595 Float_t dphistarminabs = 1e5;
596 Float_t dphistarmin = 1e5;
597
598 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
599 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
600 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
601 Float_t dphistarabs = TMath::Abs(dphistar);
602
603 if (dphistarabs < dphistarminabs) {
604 dphistarmin = dphistar;
605 dphistarminabs = dphistarabs;
606 }
607 }
608
609 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
610 //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));
611 continue;
612 }
613 }
614 }
615 fHistHBTafter->Fill(deta,dphi);
616 }//HBT cut
617
618 // conversions
619 if(fConversionCut) {
620 if (charge1 * charge2 < 0) {
621 Double_t deta = firstEta - secondEta[j];
622 Double_t dphi = firstPhi - secondPhi[j];
623 fHistConversionbefore->Fill(deta,dphi);
624
625 Float_t m0 = 0.510e-3;
626 Float_t tantheta1 = 1e10;
627
628 // phi in rad
c8ad9635 629 //Float_t phi1rad = firstPhi*TMath::DegToRad();
630 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
631 Float_t phi1rad = firstPhi;
632 Float_t phi2rad = secondPhi[j];
73609a48 633
634 if (firstEta < -1e-10 || firstEta > 1e-10)
635 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
636
637 Float_t tantheta2 = 1e10;
638 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
639 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
640
641 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
642 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
643
644 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
645
646 if (masssqu < 0.04*0.04){
647 //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));
648 continue;
649 }
650 fHistConversionafter->Fill(deta,dphi);
651 }
652 }//conversion cut
9afe3098 653
35aff0f3 654 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
655 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
656 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
657 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
8795e993 658 else {
659 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
660 continue;
661 }
6acdbcb2 662 }//end of 2nd particle loop
663 }//end of 1st particle loop
0879e280 664}
665
0879e280 666//____________________________________________________________________//
9afe3098 667TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
668 Int_t iVariablePair,
9afe3098 669 Double_t psiMin,
6acdbcb2 670 Double_t psiMax,
671 Double_t ptTriggerMin,
672 Double_t ptTriggerMax,
673 Double_t ptAssociatedMin,
674 Double_t ptAssociatedMax) {
9afe3098 675 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 676 //(private members) of the AliBalancePsi class.
621329de 677 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
678 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 679
680 // security checks
b0d53e32 681 if(psiMin > psiMax-0.00001){
682 AliError("psiMax <= psiMin");
683 return NULL;
684 }
685 if(ptTriggerMin > ptTriggerMax-0.00001){
686 AliError("ptTriggerMax <= ptTriggerMin");
687 return NULL;
688 }
689 if(ptAssociatedMin > ptAssociatedMax-0.00001){
690 AliError("ptAssociatedMax <= ptAssociatedMin");
691 return NULL;
692 }
622473b9 693
9afe3098 694 // Psi_2
622473b9 695 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
696 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
697 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
698 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
699 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
700 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
9afe3098 701
6acdbcb2 702 // pt trigger
703 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 704 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
705 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
706 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
707 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
708 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
709 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 710 }
711
712 // pt associated
713 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 714 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
715 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
716 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
717 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 718 }
719
a38fd7f3 720 //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));
721
9afe3098 722 // Project into the wanted space (1st: analysis step, 2nd: axis)
723 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
724 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
725 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
726 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
727 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
728 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
729
730 TH1D *gHistBalanceFunctionHistogram = 0x0;
731 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
732 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
733 gHistBalanceFunctionHistogram->Reset();
734
735 switch(iVariablePair) {
621329de 736 case 1:
c8ad9635 737 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
738 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 739 break;
621329de 740 case 2:
c8ad9635 741 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
742 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 743 break;
9afe3098 744 default:
745 break;
746 }
0879e280 747
0879e280 748 hTemp1->Sumw2();
749 hTemp2->Sumw2();
750 hTemp3->Sumw2();
751 hTemp4->Sumw2();
a38fd7f3 752 hTemp1->Add(hTemp3,-1.);
0879e280 753 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 754 hTemp2->Add(hTemp4,-1.);
0879e280 755 hTemp2->Scale(1./hTemp6->GetEntries());
756 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 757 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 758
759 //normalize to bin width
760 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
0879e280 761 }
762
0879e280 763 return gHistBalanceFunctionHistogram;
764}
a38fd7f3 765
f0c5040c 766//____________________________________________________________________//
767TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
768 Int_t iVariablePair,
769 Double_t psiMin,
770 Double_t psiMax,
771 Double_t ptTriggerMin,
772 Double_t ptTriggerMax,
773 Double_t ptAssociatedMin,
774 Double_t ptAssociatedMax,
775 AliBalancePsi *bfMix) {
776 //Returns the BF histogram, extracted from the 6 AliTHn objects
777 //after dividing each correlation function by the Event Mixing one
778 //(private members) of the AliBalancePsi class.
779 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
780 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 781
782 // security checks
b0d53e32 783 if(psiMin > psiMax-0.00001){
784 AliError("psiMax <= psiMin");
785 return NULL;
786 }
787 if(ptTriggerMin > ptTriggerMax-0.00001){
788 AliError("ptTriggerMax <= ptTriggerMin");
789 return NULL;
790 }
791 if(ptAssociatedMin > ptAssociatedMax-0.00001){
792 AliError("ptAssociatedMax <= ptAssociatedMin");
793 return NULL;
794 }
622473b9 795
f0c5040c 796 // Psi_2
622473b9 797 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
798 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
799 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
800 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
801 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
802 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 803
804 // pt trigger
805 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 806 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
807 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
808 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
809 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
810 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
811 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 812 }
813
814 // pt associated
815 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 816 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
817 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
818 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
819 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 820 }
821
822 // ============================================================================================
823 // the same for event mixing
824 AliTHn *fHistPMix = bfMix->GetHistNp();
825 AliTHn *fHistNMix = bfMix->GetHistNn();
826 AliTHn *fHistPNMix = bfMix->GetHistNpn();
827 AliTHn *fHistNPMix = bfMix->GetHistNnp();
828 AliTHn *fHistPPMix = bfMix->GetHistNpp();
829 AliTHn *fHistNNMix = bfMix->GetHistNnn();
830
831 // Psi_2
622473b9 832 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
833 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
834 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
835 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
836 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
837 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 838
839 // pt trigger
840 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 841 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
842 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
843 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
844 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
845 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
846 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 847 }
848
849 // pt associated
850 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 851 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
852 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
853 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
854 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 855 }
856 // ============================================================================================
857
858 //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));
859
860 // Project into the wanted space (1st: analysis step, 2nd: axis)
861 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
862 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
863 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
864 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
865 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
866 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
867
868 // ============================================================================================
869 // the same for event mixing
f0fb4ac1 870 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
871 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
872 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
873 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
874 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
875 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
f0c5040c 876 // ============================================================================================
877
878 TH1D *gHistBalanceFunctionHistogram = 0x0;
879 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
880 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
881 gHistBalanceFunctionHistogram->Reset();
882
883 switch(iVariablePair) {
884 case 1:
885 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
886 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
887 break;
888 case 2:
889 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
890 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
891 break;
892 default:
893 break;
894 }
895
896 hTemp1->Sumw2();
897 hTemp2->Sumw2();
898 hTemp3->Sumw2();
899 hTemp4->Sumw2();
900 hTemp1Mix->Sumw2();
901 hTemp2Mix->Sumw2();
902 hTemp3Mix->Sumw2();
903 hTemp4Mix->Sumw2();
904
905 hTemp1->Scale(1./hTemp5->GetEntries());
906 hTemp3->Scale(1./hTemp5->GetEntries());
907 hTemp2->Scale(1./hTemp6->GetEntries());
908 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 909 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
910 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
911 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
912 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 913
914 hTemp1->Divide(hTemp1Mix);
915 hTemp2->Divide(hTemp2Mix);
916 hTemp3->Divide(hTemp3Mix);
917 hTemp4->Divide(hTemp4Mix);
918
919 hTemp1->Add(hTemp3,-1.);
920 hTemp2->Add(hTemp4,-1.);
921
922 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
923 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 924
925 //normalize to bin width
926 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 927 }
928
929 return gHistBalanceFunctionHistogram;
930}
931
a38fd7f3 932//____________________________________________________________________//
621329de 933TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 934 Double_t psiMax,
935 Double_t ptTriggerMin,
936 Double_t ptTriggerMax,
937 Double_t ptAssociatedMin,
938 Double_t ptAssociatedMax) {
621329de 939 //Returns the BF histogram in Delta eta vs Delta phi,
940 //extracted from the 6 AliTHn objects
941 //(private members) of the AliBalancePsi class.
942 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
943 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 944
945 // security checks
b0d53e32 946 if(psiMin > psiMax-0.00001){
947 AliError("psiMax <= psiMin");
948 return NULL;
949 }
950 if(ptTriggerMin > ptTriggerMax-0.00001){
951 AliError("ptTriggerMax <= ptTriggerMin");
952 return NULL;
953 }
954 if(ptAssociatedMin > ptAssociatedMax-0.00001){
955 AliError("ptAssociatedMax <= ptAssociatedMin");
956 return NULL;
957 }
622473b9 958
621329de 959 TString histName = "gHistBalanceFunctionHistogram2D";
960
961 // Psi_2
622473b9 962 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
963 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
964 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
965 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
966 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
967 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 968
6acdbcb2 969 // pt trigger
970 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 971 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
972 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
973 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
974 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
975 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
976 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 977 }
978
979 // pt associated
980 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 981 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
982 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
983 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
984 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 985 }
986
987 //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 988
989 // Project into the wanted space (1st: analysis step, 2nd: axis)
990 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
991 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
992 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
993 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
994 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
995 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
996
997 TH2D *gHistBalanceFunctionHistogram = 0x0;
998 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
999 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1000 gHistBalanceFunctionHistogram->Reset();
c8ad9635 1001 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1002 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1003 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 1004
1005 hTemp1->Sumw2();
1006 hTemp2->Sumw2();
1007 hTemp3->Sumw2();
1008 hTemp4->Sumw2();
1009 hTemp1->Add(hTemp3,-1.);
1010 hTemp1->Scale(1./hTemp5->GetEntries());
1011 hTemp2->Add(hTemp4,-1.);
1012 hTemp2->Scale(1./hTemp6->GetEntries());
1013 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1014 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1015
1016 //normalize to bin width
1017 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
621329de 1018 }
1019
1020 return gHistBalanceFunctionHistogram;
1021}
1022
f0c5040c 1023//____________________________________________________________________//
1024TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1025 Double_t psiMax,
1026 Double_t ptTriggerMin,
1027 Double_t ptTriggerMax,
1028 Double_t ptAssociatedMin,
1029 Double_t ptAssociatedMax,
1030 AliBalancePsi *bfMix) {
1031 //Returns the BF histogram in Delta eta vs Delta phi,
1032 //after dividing each correlation function by the Event Mixing one
1033 //extracted from the 6 AliTHn objects
1034 //(private members) of the AliBalancePsi class.
1035 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1036 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1037 TString histName = "gHistBalanceFunctionHistogram2D";
1038
1039 if(!bfMix){
1040 AliError("balance function object for event mixing not available");
1041 return NULL;
1042 }
1043
622473b9 1044 // security checks
b0d53e32 1045 if(psiMin > psiMax-0.00001){
1046 AliError("psiMax <= psiMin");
1047 return NULL;
1048 }
1049 if(ptTriggerMin > ptTriggerMax-0.00001){
1050 AliError("ptTriggerMax <= ptTriggerMin");
1051 return NULL;
1052 }
1053 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1054 AliError("ptAssociatedMax <= ptAssociatedMin");
1055 return NULL;
1056 }
622473b9 1057
f0c5040c 1058 // Psi_2
622473b9 1059 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1060 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1061 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1062 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1063 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1064 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1065
1066 // pt trigger
1067 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1068 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1069 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1070 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1071 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1072 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1073 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1074 }
1075
1076 // pt associated
1077 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1078 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1079 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1080 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1081 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1082 }
1083
1084
1085 // ============================================================================================
1086 // the same for event mixing
1087 AliTHn *fHistPMix = bfMix->GetHistNp();
1088 AliTHn *fHistNMix = bfMix->GetHistNn();
1089 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1090 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1091 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1092 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1093
1094 // Psi_2
622473b9 1095 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1096 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1097 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1098 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1099 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1100 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1101
1102 // pt trigger
1103 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1104 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1105 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1106 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1107 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1108 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1109 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1110 }
1111
1112 // pt associated
1113 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1114 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1115 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1116 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1117 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1118 }
1119 // ============================================================================================
1120
1121
1122 //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)));
1123
1124 // Project into the wanted space (1st: analysis step, 2nd: axis)
1125 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1126 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1127 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1128 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1129 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1130 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1131
1132 // ============================================================================================
1133 // the same for event mixing
1134 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1135 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1136 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1137 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
f0fb4ac1 1138 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1139 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
f0c5040c 1140 // ============================================================================================
1141
1142 TH2D *gHistBalanceFunctionHistogram = 0x0;
1143 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1144 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1145 gHistBalanceFunctionHistogram->Reset();
1146 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1147 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1148 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1149
1150 hTemp1->Sumw2();
1151 hTemp2->Sumw2();
1152 hTemp3->Sumw2();
1153 hTemp4->Sumw2();
1154 hTemp1Mix->Sumw2();
1155 hTemp2Mix->Sumw2();
1156 hTemp3Mix->Sumw2();
1157 hTemp4Mix->Sumw2();
1158
1159 hTemp1->Scale(1./hTemp5->GetEntries());
1160 hTemp3->Scale(1./hTemp5->GetEntries());
1161 hTemp2->Scale(1./hTemp6->GetEntries());
1162 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 1163 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1164 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1165 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1166 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 1167
1168 hTemp1->Divide(hTemp1Mix);
1169 hTemp2->Divide(hTemp2Mix);
1170 hTemp3->Divide(hTemp3Mix);
1171 hTemp4->Divide(hTemp4Mix);
1172
1173 hTemp1->Add(hTemp3,-1.);
1174 hTemp2->Add(hTemp4,-1.);
1175
1176 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1177 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1178
1179 //normalize to bin width
1180 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 1181 }
1182
1183 return gHistBalanceFunctionHistogram;
1184}
1185
f2e8af26 1186//____________________________________________________________________//
1187TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1188 Double_t psiMin,
1189 Double_t psiMax,
1190 Double_t ptTriggerMin,
1191 Double_t ptTriggerMax,
1192 Double_t ptAssociatedMin,
1193 Double_t ptAssociatedMax,
1194 AliBalancePsi *bfMix) {
1195 //Returns the BF histogram in Delta eta OR Delta phi,
1196 //after dividing each correlation function by the Event Mixing one
1197 // (But the division is done here in 2D, this was basically done to check the Integral)
1198 //extracted from the 6 AliTHn objects
1199 //(private members) of the AliBalancePsi class.
1200 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1201 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1202 TString histName = "gHistBalanceFunctionHistogram1D";
1203
1204 if(!bfMix){
1205 AliError("balance function object for event mixing not available");
1206 return NULL;
1207 }
1208
622473b9 1209 // security checks
b0d53e32 1210 if(psiMin > psiMax-0.00001){
1211 AliError("psiMax <= psiMin");
1212 return NULL;
1213 }
1214 if(ptTriggerMin > ptTriggerMax-0.00001){
1215 AliError("ptTriggerMax <= ptTriggerMin");
1216 return NULL;
1217 }
1218 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1219 AliError("ptAssociatedMax <= ptAssociatedMin");
1220 return NULL;
1221 }
622473b9 1222
f2e8af26 1223 // Psi_2
622473b9 1224 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1225 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1226 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1227 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1228 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1229 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1230
1231 // pt trigger
1232 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1233 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1234 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1235 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1236 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1237 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1238 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1239 }
1240
1241 // pt associated
1242 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1243 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1244 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1245 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1246 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1247 }
1248
1249
1250 // ============================================================================================
1251 // the same for event mixing
1252 AliTHn *fHistPMix = bfMix->GetHistNp();
1253 AliTHn *fHistNMix = bfMix->GetHistNn();
1254 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1255 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1256 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1257 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1258
1259 // Psi_2
622473b9 1260 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1261 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1262 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1263 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1264 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1265 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1266
1267 // pt trigger
1268 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1269 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1270 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1271 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1272 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1273 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1274 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1275 }
1276
1277 // pt associated
1278 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1279 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1280 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1281 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1282 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1283 }
1284 // ============================================================================================
1285
1286
1287 //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)));
1288
1289 // Project into the wanted space (1st: analysis step, 2nd: axis)
1290 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1291 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1292 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1293 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1294 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1295 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1296
1297 // ============================================================================================
1298 // the same for event mixing
1299 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1300 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1301 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1302 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1303 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1304 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1305 // ============================================================================================
1306
1307 TH1D *gHistBalanceFunctionHistogram = 0x0;
1308 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1309
1310 hTemp1->Sumw2();
1311 hTemp2->Sumw2();
1312 hTemp3->Sumw2();
1313 hTemp4->Sumw2();
1314 hTemp1Mix->Sumw2();
1315 hTemp2Mix->Sumw2();
1316 hTemp3Mix->Sumw2();
1317 hTemp4Mix->Sumw2();
1318
1319 hTemp1->Scale(1./hTemp5->GetEntries());
1320 hTemp3->Scale(1./hTemp5->GetEntries());
1321 hTemp2->Scale(1./hTemp6->GetEntries());
1322 hTemp4->Scale(1./hTemp6->GetEntries());
1323
1324 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1325 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1326 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1327 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
1328
1329 hTemp1->Divide(hTemp1Mix);
1330 hTemp2->Divide(hTemp2Mix);
1331 hTemp3->Divide(hTemp3Mix);
1332 hTemp4->Divide(hTemp4Mix);
1333
1334 // now only project on one axis
1335 TH1D *h1DTemp1 = NULL;
1336 TH1D *h1DTemp2 = NULL;
1337 TH1D *h1DTemp3 = NULL;
1338 TH1D *h1DTemp4 = NULL;
1339 if(!bPhi){
1340 h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
1341 h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
1342 h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
1343 h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
1344 }
1345 else{
1346 h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
1347 h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
1348 h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
1349 h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
1350 }
1351
1352 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1353 gHistBalanceFunctionHistogram->Reset();
1354 if(!bPhi){
1355 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1356 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1357 }
1358 else{
1359 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1360 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1361 }
1362
1363 h1DTemp1->Add(h1DTemp3,-1.);
1364 h1DTemp2->Add(h1DTemp4,-1.);
1365
1366 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1367 gHistBalanceFunctionHistogram->Scale(0.5);
1368 }
1369
1370 return gHistBalanceFunctionHistogram;
1371}
1372
f0c5040c 1373
621329de 1374//____________________________________________________________________//
1375TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1376 Double_t psiMax,
1377 Double_t ptTriggerMin,
1378 Double_t ptTriggerMax,
1379 Double_t ptAssociatedMin,
47e040b5 1380 Double_t ptAssociatedMax) {
a38fd7f3 1381 //Returns the 2D correlation function for (+-) pairs
622473b9 1382
1383 // security checks
b0d53e32 1384 if(psiMin > psiMax-0.00001){
1385 AliError("psiMax <= psiMin");
1386 return NULL;
1387 }
1388 if(ptTriggerMin > ptTriggerMax-0.00001){
1389 AliError("ptTriggerMax <= ptTriggerMin");
1390 return NULL;
1391 }
1392 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1393 AliError("ptAssociatedMax <= ptAssociatedMin");
1394 return NULL;
1395 }
622473b9 1396
621329de 1397 // Psi_2: axis 0
622473b9 1398 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1399 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1400
6acdbcb2 1401 // pt trigger
1402 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1403 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1404 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1405 }
621329de 1406
6acdbcb2 1407 // pt associated
1408 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1409 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1410
1411 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1412 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1413
1414 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1415 //TCanvas *c1 = new TCanvas("c1","");
1416 //c1->cd();
1417 //if(!gHistTest){
1418 //AliError("Projection of fHistP = NULL");
1419 //return gHistTest;
1420 //}
1421 //else{
1422 //gHistTest->DrawCopy("colz");
1423 //}
1424
1425 //0:step, 1: Delta eta, 2: Delta phi
621329de 1426 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1427 if(!gHist){
1428 AliError("Projection of fHistPN = NULL");
1429 return gHist;
1430 }
1431
6acdbcb2 1432 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1433 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1434 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1435
6acdbcb2 1436 //TCanvas *c2 = new TCanvas("c2","");
1437 //c2->cd();
1438 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1439
47e040b5 1440 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1441 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1442
1443 //normalize to bin width
1444 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1445
1446 return gHist;
1447}
1448
1449//____________________________________________________________________//
621329de 1450TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1451 Double_t psiMax,
1452 Double_t ptTriggerMin,
1453 Double_t ptTriggerMax,
1454 Double_t ptAssociatedMin,
47e040b5 1455 Double_t ptAssociatedMax) {
a38fd7f3 1456 //Returns the 2D correlation function for (+-) pairs
622473b9 1457
1458 // security checks
b0d53e32 1459 if(psiMin > psiMax-0.00001){
1460 AliError("psiMax <= psiMin");
1461 return NULL;
1462 }
1463 if(ptTriggerMin > ptTriggerMax-0.00001){
1464 AliError("ptTriggerMax <= ptTriggerMin");
1465 return NULL;
1466 }
1467 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1468 AliError("ptAssociatedMax <= ptAssociatedMin");
1469 return NULL;
1470 }
622473b9 1471
621329de 1472 // Psi_2: axis 0
622473b9 1473 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1474 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
a38fd7f3 1475
6acdbcb2 1476 // pt trigger
1477 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1478 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1479 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1480 }
1481
1482 // pt associated
1483 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1484 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1485
1486 //0:step, 1: Delta eta, 2: Delta phi
621329de 1487 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1488 if(!gHist){
1489 AliError("Projection of fHistPN = NULL");
1490 return gHist;
1491 }
1492
a38fd7f3 1493 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1494 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 1495 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1496 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1497
1498 //normalize to bin width
1499 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1500
1501 return gHist;
1502}
1503
1504//____________________________________________________________________//
621329de 1505TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1506 Double_t psiMax,
1507 Double_t ptTriggerMin,
1508 Double_t ptTriggerMax,
1509 Double_t ptAssociatedMin,
47e040b5 1510 Double_t ptAssociatedMax) {
a38fd7f3 1511 //Returns the 2D correlation function for (+-) pairs
622473b9 1512
1513 // security checks
b0d53e32 1514 if(psiMin > psiMax-0.00001){
1515 AliError("psiMax <= psiMin");
1516 return NULL;
1517 }
1518 if(ptTriggerMin > ptTriggerMax-0.00001){
1519 AliError("ptTriggerMax <= ptTriggerMin");
1520 return NULL;
1521 }
1522 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1523 AliError("ptAssociatedMax <= ptAssociatedMin");
1524 return NULL;
1525 }
622473b9 1526
621329de 1527 // Psi_2: axis 0
622473b9 1528 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1529 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1530
1531 // pt trigger
1532 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1533 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1534 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1535 }
1536
1537 // pt associated
1538 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1539 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1540
6acdbcb2 1541 //0:step, 1: Delta eta, 2: Delta phi
621329de 1542 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1543 if(!gHist){
1544 AliError("Projection of fHistPN = NULL");
1545 return gHist;
1546 }
1547
a38fd7f3 1548 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1549 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 1550 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1551 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1552
1553 //normalize to bin width
1554 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1555
1556 return gHist;
1557}
1558
1559//____________________________________________________________________//
621329de 1560TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1561 Double_t psiMax,
1562 Double_t ptTriggerMin,
1563 Double_t ptTriggerMax,
1564 Double_t ptAssociatedMin,
47e040b5 1565 Double_t ptAssociatedMax) {
a38fd7f3 1566 //Returns the 2D correlation function for (+-) pairs
622473b9 1567
1568 // security checks
1569 if(psiMin > psiMax-0.00001){
1570 AliError("psiMax <= psiMin");
1571 return NULL;
1572 }
1573 if(ptTriggerMin > ptTriggerMax-0.00001){
1574 AliError("ptTriggerMax <= ptTriggerMin");
1575 return NULL;
1576 }
1577 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1578 AliError("ptAssociatedMax <= ptAssociatedMin");
1579 return NULL;
1580 }
1581
621329de 1582 // Psi_2: axis 0
622473b9 1583 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1584 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1585
1586 // pt trigger
1587 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1588 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1589 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1590 }
1591
1592 // pt associated
1593 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1594 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1595
6acdbcb2 1596 //0:step, 1: Delta eta, 2: Delta phi
621329de 1597 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 1598 if(!gHist){
1599 AliError("Projection of fHistPN = NULL");
1600 return gHist;
1601 }
1602
a38fd7f3 1603 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1604 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 1605 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1606 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1607
1608 //normalize to bin width
1609 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1610
1611 return gHist;
1612}
621329de 1613
7b610f27 1614//____________________________________________________________________//
1615TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1616 Double_t psiMax,
1617 Double_t ptTriggerMin,
1618 Double_t ptTriggerMax,
1619 Double_t ptAssociatedMin,
47e040b5 1620 Double_t ptAssociatedMax) {
7b610f27 1621 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 1622
1623 // security checks
b0d53e32 1624 if(psiMin > psiMax-0.00001){
1625 AliError("psiMax <= psiMin");
1626 return NULL;
1627 }
1628 if(ptTriggerMin > ptTriggerMax-0.00001){
1629 AliError("ptTriggerMax <= ptTriggerMin");
1630 return NULL;
1631 }
1632 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1633 AliError("ptAssociatedMax <= ptAssociatedMin");
1634 return NULL;
1635 }
622473b9 1636
7b610f27 1637 // Psi_2: axis 0
622473b9 1638 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1639 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1640 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1641 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1642 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1643 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
7b610f27 1644
1645 // pt trigger
1646 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1647 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1648 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1649 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1650 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1651 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1652 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 1653 }
1654
1655 // pt associated
1656 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1657 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1658 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1659 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1660 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 1661
1662 //0:step, 1: Delta eta, 2: Delta phi
1663 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1664 if(!gHistNN){
1665 AliError("Projection of fHistNN = NULL");
1666 return gHistNN;
1667 }
1668 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1669 if(!gHistPP){
1670 AliError("Projection of fHistPP = NULL");
1671 return gHistPP;
1672 }
1673 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1674 if(!gHistNP){
1675 AliError("Projection of fHistNP = NULL");
1676 return gHistNP;
1677 }
1678 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1679 if(!gHistPN){
1680 AliError("Projection of fHistPN = NULL");
1681 return gHistPN;
1682 }
1683
1684 // sum all 2 particle histograms
1685 gHistNN->Add(gHistPP);
1686 gHistNN->Add(gHistNP);
1687 gHistNN->Add(gHistPN);
1688
47e040b5 1689 // divide by sum of + and - triggers
1690 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 1691 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 1692
1693 //normalize to bin width
1694 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 1695
1696 return gHistNN;
1697}
1698
f2e8af26 1699//____________________________________________________________________//
1700
1701Bool_t AliBalancePsi::GetMomentsAnalytical(TH1D* gHist,
1702 Double_t &mean, Double_t &meanError,
1703 Double_t &sigma, Double_t &sigmaError,
1704 Double_t &skewness, Double_t &skewnessError,
1705 Double_t &kurtosis, Double_t &kurtosisError) {
1706 //
1707 // helper method to calculate the moments and errors of a TH1D anlytically
1708 //
1709
1710 Bool_t success = kFALSE;
1711 mean = -1.;
1712 meanError = -1.;
1713 sigma = -1.;
1714 sigmaError = -1.;
1715 skewness = -1.;
1716 skewnessError = -1.;
1717 kurtosis = -1.;
1718 kurtosisError = -1.;
1719
1720 if(gHist){
1721
1722 // ----------------------------------------------------------------------
1723 // basic parameters of histogram
1724
1725 Int_t fNumberOfBins = gHist->GetNbinsX();
1726 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
1727
1728
1729 // ----------------------------------------------------------------------
1730 // first calculate the mean
1731
1732 Double_t fWeightedAverage = 0.;
f0754617 1733 Double_t fNormalization = 0.;
f2e8af26 1734
1735 for(Int_t i = 1; i <= fNumberOfBins; i++) {
1736
1737 fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
1738 fNormalization += gHist->GetBinContent(i);
f2e8af26 1739 }
1740
1741 mean = fWeightedAverage / fNormalization;
1742
f2e8af26 1743 // ----------------------------------------------------------------------
1744 // then calculate the higher moments
1745
1746 Double_t fDelta = 0.;
1747 Double_t fDelta2 = 0.;
1748 Double_t fDelta3 = 0.;
1749 Double_t fDelta4 = 0.;
1750
1751 for(Int_t i = 1; i <= fNumberOfBins; i++) {
f0754617 1752 fDelta += gHist->GetBinContent(i) * (gHist->GetBinCenter(i) - mean);
1753 fDelta2 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
1754 fDelta3 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
1755 fDelta4 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
f2e8af26 1756 }
1757
f0754617 1758 sigma = TMath::Sqrt(fDelta2 / fNormalization);
f2e8af26 1759 skewness = fDelta3 / fNormalization / TMath::Power(sigma,3);
1760 kurtosis = fDelta4 / fNormalization / TMath::Power(sigma,4) - 3;
1761
f0754617 1762 // ----------------------------------------------------------------------
1763 // then calculate the higher moment errors?
1764 cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
1765 meanError = sigma / TMath::Sqrt(fNormalization); //sigma / TMath::Sqrt(gHist->GetEffectiveEntries()) would give the "ROOT" menaError
1766
f2e8af26 1767 success = kTRUE;
1768 }
1769
1770
1771 return success;
1772}
1773
1774
73609a48 1775//____________________________________________________________________//
1776Float_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) {
1777 //
1778 // calculates dphistar
1779 //
1780 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1781
1782 static const Double_t kPi = TMath::Pi();
1783
1784 // circularity
1785// if (dphistar > 2 * kPi)
1786// dphistar -= 2 * kPi;
1787// if (dphistar < -2 * kPi)
1788// dphistar += 2 * kPi;
1789
1790 if (dphistar > kPi)
1791 dphistar = kPi * 2 - dphistar;
1792 if (dphistar < -kPi)
1793 dphistar = -kPi * 2 - dphistar;
1794 if (dphistar > kPi) // might look funny but is needed
1795 dphistar = kPi * 2 - dphistar;
1796
1797 return dphistar;
1798}
1799
2922bbe3 1800
1801