]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliBalancePsi.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updated
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalancePsi.cxx
CommitLineData
0879e280 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14/* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
15
16//-----------------------------------------------------------------
17// Balance Function class
18// This is the class to deal with the Balance Function wrt Psi analysis
19// Origin: Panos Christakoglou, Nikhef, Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21
22
23//ROOT
24#include <Riostream.h>
621329de 25#include <TCanvas.h>
0879e280 26#include <TMath.h>
27#include <TAxis.h>
28#include <TH2D.h>
29#include <TH3D.h>
30#include <TLorentzVector.h>
31#include <TObjArray.h>
32#include <TGraphErrors.h>
33#include <TString.h>
34
35#include "AliVParticle.h"
36#include "AliMCParticle.h"
37#include "AliESDtrack.h"
38#include "AliAODTrack.h"
9afe3098 39#include "AliTHn.h"
35aff0f3 40#include "AliAnalysisTaskTriggeredBF.h"
0879e280 41
42#include "AliBalancePsi.h"
43
44ClassImp(AliBalancePsi)
45
46//____________________________________________________________________//
47AliBalancePsi::AliBalancePsi() :
48 TObject(),
9fd4b54e 49 fShuffle(kFALSE),
0879e280 50 fAnalysisLevel("ESD"),
51 fAnalyzedEvents(0) ,
52 fCentralityId(0) ,
53 fCentStart(0.),
54 fCentStop(0.),
9afe3098 55 fHistP(0),
56 fHistN(0),
57 fHistPN(0),
58 fHistNP(0),
59 fHistPP(0),
60 fHistNN(0),
73609a48 61 fHistHBTbefore(0),
62 fHistHBTafter(0),
63 fHistConversionbefore(0),
64 fHistConversionafter(0),
c8ad9635 65 fHistPsiMinusPhi(0),
f0754617 66 fHistResonancesBefore(0),
67 fHistResonancesRho(0),
68 fHistResonancesK0(0),
69 fHistResonancesLambda(0),
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) {
749de081 535 if (charge1 * charge2 < 0) {
536
537 //rho0
538 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
539 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
540 vectorMother = vectorDaughter[0] + vectorDaughter[1];
541 fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
542 if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0)
543 continue;
544 fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
545
546 //K0s
547 if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s)
548 continue;
549 fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
550
551
552 //Lambda
553 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
554 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass());
555 vectorMother = vectorDaughter[0] + vectorDaughter[1];
556 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
557 continue;
558
559 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass());
560 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
561 vectorMother = vectorDaughter[0] + vectorDaughter[1];
562 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
563 continue;
564 fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
565
566 }//unlike-sign only
2922bbe3 567 }//resonance cut
73609a48 568
569 // HBT like cut
da8b2d30 570 if(fHBTCut){ // VERSION 3 (all pairs)
571 //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
73609a48 572 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
573 // continue;
574
575 Double_t deta = firstEta - secondEta[j];
576 Double_t dphi = firstPhi - secondPhi[j];
577 // VERSION 2 (Taken from DPhiCorrelations)
578 // the variables & cuthave been developed by the HBT group
579 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
580 fHistHBTbefore->Fill(deta,dphi);
581
582 // optimization
583 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
584 {
585 // phi in rad
c8ad9635 586 //Float_t phi1rad = firstPhi*TMath::DegToRad();
587 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
588 Float_t phi1rad = firstPhi;
589 Float_t phi2rad = secondPhi[j];
73609a48 590
591 // check first boundaries to see if is worth to loop and find the minimum
592 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
593 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
594
595 const Float_t kLimit = 0.02 * 3;
596
597 Float_t dphistarminabs = 1e5;
598 Float_t dphistarmin = 1e5;
599
600 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
601 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
602 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
603 Float_t dphistarabs = TMath::Abs(dphistar);
604
605 if (dphistarabs < dphistarminabs) {
606 dphistarmin = dphistar;
607 dphistarminabs = dphistarabs;
608 }
609 }
610
611 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
612 //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));
613 continue;
614 }
615 }
616 }
617 fHistHBTafter->Fill(deta,dphi);
618 }//HBT cut
619
620 // conversions
621 if(fConversionCut) {
622 if (charge1 * charge2 < 0) {
623 Double_t deta = firstEta - secondEta[j];
624 Double_t dphi = firstPhi - secondPhi[j];
625 fHistConversionbefore->Fill(deta,dphi);
626
627 Float_t m0 = 0.510e-3;
628 Float_t tantheta1 = 1e10;
629
630 // phi in rad
c8ad9635 631 //Float_t phi1rad = firstPhi*TMath::DegToRad();
632 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
633 Float_t phi1rad = firstPhi;
634 Float_t phi2rad = secondPhi[j];
73609a48 635
636 if (firstEta < -1e-10 || firstEta > 1e-10)
637 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
638
639 Float_t tantheta2 = 1e10;
640 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
641 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
642
643 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
644 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
645
646 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
647
648 if (masssqu < 0.04*0.04){
649 //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));
650 continue;
651 }
652 fHistConversionafter->Fill(deta,dphi);
653 }
654 }//conversion cut
9afe3098 655
35aff0f3 656 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
657 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
658 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
659 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
8795e993 660 else {
661 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
662 continue;
663 }
6acdbcb2 664 }//end of 2nd particle loop
665 }//end of 1st particle loop
0879e280 666}
667
0879e280 668//____________________________________________________________________//
9afe3098 669TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
670 Int_t iVariablePair,
9afe3098 671 Double_t psiMin,
6acdbcb2 672 Double_t psiMax,
673 Double_t ptTriggerMin,
674 Double_t ptTriggerMax,
675 Double_t ptAssociatedMin,
676 Double_t ptAssociatedMax) {
9afe3098 677 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 678 //(private members) of the AliBalancePsi class.
621329de 679 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
680 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 681
682 // security checks
b0d53e32 683 if(psiMin > psiMax-0.00001){
684 AliError("psiMax <= psiMin");
685 return NULL;
686 }
687 if(ptTriggerMin > ptTriggerMax-0.00001){
688 AliError("ptTriggerMax <= ptTriggerMin");
689 return NULL;
690 }
691 if(ptAssociatedMin > ptAssociatedMax-0.00001){
692 AliError("ptAssociatedMax <= ptAssociatedMin");
693 return NULL;
694 }
622473b9 695
9afe3098 696 // Psi_2
622473b9 697 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
698 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
699 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
700 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
701 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
702 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
9afe3098 703
6acdbcb2 704 // pt trigger
705 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 706 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
707 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
708 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
709 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
710 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
711 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 712 }
713
714 // pt associated
715 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 716 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
717 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
718 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
719 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 720 }
721
a38fd7f3 722 //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));
723
9afe3098 724 // Project into the wanted space (1st: analysis step, 2nd: axis)
725 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
726 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
727 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
728 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
729 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
730 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
731
732 TH1D *gHistBalanceFunctionHistogram = 0x0;
733 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
734 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
735 gHistBalanceFunctionHistogram->Reset();
736
737 switch(iVariablePair) {
621329de 738 case 1:
c8ad9635 739 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
740 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 741 break;
621329de 742 case 2:
c8ad9635 743 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
744 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 745 break;
9afe3098 746 default:
747 break;
748 }
0879e280 749
0879e280 750 hTemp1->Sumw2();
751 hTemp2->Sumw2();
752 hTemp3->Sumw2();
753 hTemp4->Sumw2();
a38fd7f3 754 hTemp1->Add(hTemp3,-1.);
0879e280 755 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 756 hTemp2->Add(hTemp4,-1.);
0879e280 757 hTemp2->Scale(1./hTemp6->GetEntries());
758 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 759 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 760
761 //normalize to bin width
762 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
0879e280 763 }
764
0879e280 765 return gHistBalanceFunctionHistogram;
766}
a38fd7f3 767
f0c5040c 768//____________________________________________________________________//
769TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
770 Int_t iVariablePair,
771 Double_t psiMin,
772 Double_t psiMax,
773 Double_t ptTriggerMin,
774 Double_t ptTriggerMax,
775 Double_t ptAssociatedMin,
776 Double_t ptAssociatedMax,
777 AliBalancePsi *bfMix) {
778 //Returns the BF histogram, extracted from the 6 AliTHn objects
779 //after dividing each correlation function by the Event Mixing one
780 //(private members) of the AliBalancePsi class.
781 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
782 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 783
784 // security checks
b0d53e32 785 if(psiMin > psiMax-0.00001){
786 AliError("psiMax <= psiMin");
787 return NULL;
788 }
789 if(ptTriggerMin > ptTriggerMax-0.00001){
790 AliError("ptTriggerMax <= ptTriggerMin");
791 return NULL;
792 }
793 if(ptAssociatedMin > ptAssociatedMax-0.00001){
794 AliError("ptAssociatedMax <= ptAssociatedMin");
795 return NULL;
796 }
622473b9 797
f0c5040c 798 // Psi_2
622473b9 799 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
800 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
801 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
802 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
803 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
804 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 805
806 // pt trigger
807 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 808 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
809 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
810 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
811 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
812 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
813 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 814 }
815
816 // pt associated
817 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 818 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
819 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
820 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
821 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 822 }
823
824 // ============================================================================================
825 // the same for event mixing
826 AliTHn *fHistPMix = bfMix->GetHistNp();
827 AliTHn *fHistNMix = bfMix->GetHistNn();
828 AliTHn *fHistPNMix = bfMix->GetHistNpn();
829 AliTHn *fHistNPMix = bfMix->GetHistNnp();
830 AliTHn *fHistPPMix = bfMix->GetHistNpp();
831 AliTHn *fHistNNMix = bfMix->GetHistNnn();
832
833 // Psi_2
622473b9 834 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
835 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
836 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
837 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
838 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
839 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 840
841 // pt trigger
842 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 843 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
844 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
845 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
846 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
847 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
848 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 849 }
850
851 // pt associated
852 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 853 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
854 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
855 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
856 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 857 }
858 // ============================================================================================
859
860 //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));
861
862 // Project into the wanted space (1st: analysis step, 2nd: axis)
863 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
864 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
865 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
866 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
867 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
868 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
869
870 // ============================================================================================
871 // the same for event mixing
f0fb4ac1 872 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
873 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
874 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
875 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
876 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
877 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
f0c5040c 878 // ============================================================================================
879
880 TH1D *gHistBalanceFunctionHistogram = 0x0;
881 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
882 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
883 gHistBalanceFunctionHistogram->Reset();
884
885 switch(iVariablePair) {
886 case 1:
887 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
888 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
889 break;
890 case 2:
891 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
892 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
893 break;
894 default:
895 break;
896 }
897
898 hTemp1->Sumw2();
899 hTemp2->Sumw2();
900 hTemp3->Sumw2();
901 hTemp4->Sumw2();
902 hTemp1Mix->Sumw2();
903 hTemp2Mix->Sumw2();
904 hTemp3Mix->Sumw2();
905 hTemp4Mix->Sumw2();
906
907 hTemp1->Scale(1./hTemp5->GetEntries());
908 hTemp3->Scale(1./hTemp5->GetEntries());
909 hTemp2->Scale(1./hTemp6->GetEntries());
910 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 911 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
912 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
913 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
914 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 915
916 hTemp1->Divide(hTemp1Mix);
917 hTemp2->Divide(hTemp2Mix);
918 hTemp3->Divide(hTemp3Mix);
919 hTemp4->Divide(hTemp4Mix);
920
921 hTemp1->Add(hTemp3,-1.);
922 hTemp2->Add(hTemp4,-1.);
923
924 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
925 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 926
927 //normalize to bin width
928 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 929 }
930
931 return gHistBalanceFunctionHistogram;
932}
933
a38fd7f3 934//____________________________________________________________________//
621329de 935TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 936 Double_t psiMax,
937 Double_t ptTriggerMin,
938 Double_t ptTriggerMax,
939 Double_t ptAssociatedMin,
940 Double_t ptAssociatedMax) {
621329de 941 //Returns the BF histogram in Delta eta vs Delta phi,
942 //extracted from the 6 AliTHn objects
943 //(private members) of the AliBalancePsi class.
944 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
945 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 946
947 // security checks
b0d53e32 948 if(psiMin > psiMax-0.00001){
949 AliError("psiMax <= psiMin");
950 return NULL;
951 }
952 if(ptTriggerMin > ptTriggerMax-0.00001){
953 AliError("ptTriggerMax <= ptTriggerMin");
954 return NULL;
955 }
956 if(ptAssociatedMin > ptAssociatedMax-0.00001){
957 AliError("ptAssociatedMax <= ptAssociatedMin");
958 return NULL;
959 }
622473b9 960
621329de 961 TString histName = "gHistBalanceFunctionHistogram2D";
962
963 // Psi_2
622473b9 964 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
965 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
966 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
967 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
968 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
969 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 970
6acdbcb2 971 // pt trigger
972 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 973 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
974 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
975 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
976 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
977 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
978 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 979 }
980
981 // pt associated
982 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 983 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
984 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
985 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
986 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 987 }
988
989 //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 990
991 // Project into the wanted space (1st: analysis step, 2nd: axis)
992 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
993 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
994 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
995 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
996 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
997 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
998
999 TH2D *gHistBalanceFunctionHistogram = 0x0;
1000 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1001 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1002 gHistBalanceFunctionHistogram->Reset();
c8ad9635 1003 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1004 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1005 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 1006
1007 hTemp1->Sumw2();
1008 hTemp2->Sumw2();
1009 hTemp3->Sumw2();
1010 hTemp4->Sumw2();
1011 hTemp1->Add(hTemp3,-1.);
1012 hTemp1->Scale(1./hTemp5->GetEntries());
1013 hTemp2->Add(hTemp4,-1.);
1014 hTemp2->Scale(1./hTemp6->GetEntries());
1015 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1016 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1017
1018 //normalize to bin width
1019 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
621329de 1020 }
1021
1022 return gHistBalanceFunctionHistogram;
1023}
1024
f0c5040c 1025//____________________________________________________________________//
1026TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1027 Double_t psiMax,
1028 Double_t ptTriggerMin,
1029 Double_t ptTriggerMax,
1030 Double_t ptAssociatedMin,
1031 Double_t ptAssociatedMax,
1032 AliBalancePsi *bfMix) {
1033 //Returns the BF histogram in Delta eta vs Delta phi,
1034 //after dividing each correlation function by the Event Mixing one
1035 //extracted from the 6 AliTHn objects
1036 //(private members) of the AliBalancePsi class.
1037 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1038 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1039 TString histName = "gHistBalanceFunctionHistogram2D";
1040
1041 if(!bfMix){
1042 AliError("balance function object for event mixing not available");
1043 return NULL;
1044 }
1045
622473b9 1046 // security checks
b0d53e32 1047 if(psiMin > psiMax-0.00001){
1048 AliError("psiMax <= psiMin");
1049 return NULL;
1050 }
1051 if(ptTriggerMin > ptTriggerMax-0.00001){
1052 AliError("ptTriggerMax <= ptTriggerMin");
1053 return NULL;
1054 }
1055 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1056 AliError("ptAssociatedMax <= ptAssociatedMin");
1057 return NULL;
1058 }
622473b9 1059
f0c5040c 1060 // Psi_2
622473b9 1061 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1062 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1063 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1064 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1065 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1066 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1067
1068 // pt trigger
1069 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1070 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1071 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1072 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1073 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1074 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1075 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1076 }
1077
1078 // pt associated
1079 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1080 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1081 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1082 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1083 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1084 }
1085
1086
1087 // ============================================================================================
1088 // the same for event mixing
1089 AliTHn *fHistPMix = bfMix->GetHistNp();
1090 AliTHn *fHistNMix = bfMix->GetHistNn();
1091 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1092 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1093 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1094 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1095
1096 // Psi_2
622473b9 1097 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1098 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1099 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1100 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1101 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1102 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1103
1104 // pt trigger
1105 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1106 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1107 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1108 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1109 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1110 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1111 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1112 }
1113
1114 // pt associated
1115 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1116 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1117 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1118 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1119 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1120 }
1121 // ============================================================================================
1122
1123
1124 //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
1125
1126 // Project into the wanted space (1st: analysis step, 2nd: axis)
1127 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1128 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1129 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1130 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1131 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1132 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1133
1134 // ============================================================================================
1135 // the same for event mixing
1136 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1137 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1138 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1139 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
f0fb4ac1 1140 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1141 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
f0c5040c 1142 // ============================================================================================
1143
1144 TH2D *gHistBalanceFunctionHistogram = 0x0;
1145 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1146 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1147 gHistBalanceFunctionHistogram->Reset();
1148 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1149 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1150 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1151
1152 hTemp1->Sumw2();
1153 hTemp2->Sumw2();
1154 hTemp3->Sumw2();
1155 hTemp4->Sumw2();
1156 hTemp1Mix->Sumw2();
1157 hTemp2Mix->Sumw2();
1158 hTemp3Mix->Sumw2();
1159 hTemp4Mix->Sumw2();
1160
1161 hTemp1->Scale(1./hTemp5->GetEntries());
1162 hTemp3->Scale(1./hTemp5->GetEntries());
1163 hTemp2->Scale(1./hTemp6->GetEntries());
1164 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 1165 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1166 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1167 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1168 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 1169
1170 hTemp1->Divide(hTemp1Mix);
1171 hTemp2->Divide(hTemp2Mix);
1172 hTemp3->Divide(hTemp3Mix);
1173 hTemp4->Divide(hTemp4Mix);
1174
1175 hTemp1->Add(hTemp3,-1.);
1176 hTemp2->Add(hTemp4,-1.);
1177
1178 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1179 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1180
1181 //normalize to bin width
1182 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 1183 }
1184
1185 return gHistBalanceFunctionHistogram;
1186}
1187
f2e8af26 1188//____________________________________________________________________//
1189TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1190 Double_t psiMin,
1191 Double_t psiMax,
1192 Double_t ptTriggerMin,
1193 Double_t ptTriggerMax,
1194 Double_t ptAssociatedMin,
1195 Double_t ptAssociatedMax,
1196 AliBalancePsi *bfMix) {
1197 //Returns the BF histogram in Delta eta OR Delta phi,
1198 //after dividing each correlation function by the Event Mixing one
1199 // (But the division is done here in 2D, this was basically done to check the Integral)
1200 //extracted from the 6 AliTHn objects
1201 //(private members) of the AliBalancePsi class.
1202 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1203 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1204 TString histName = "gHistBalanceFunctionHistogram1D";
1205
1206 if(!bfMix){
1207 AliError("balance function object for event mixing not available");
1208 return NULL;
1209 }
1210
622473b9 1211 // security checks
b0d53e32 1212 if(psiMin > psiMax-0.00001){
1213 AliError("psiMax <= psiMin");
1214 return NULL;
1215 }
1216 if(ptTriggerMin > ptTriggerMax-0.00001){
1217 AliError("ptTriggerMax <= ptTriggerMin");
1218 return NULL;
1219 }
1220 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1221 AliError("ptAssociatedMax <= ptAssociatedMin");
1222 return NULL;
1223 }
622473b9 1224
f2e8af26 1225 // Psi_2
622473b9 1226 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1227 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1228 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1229 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1230 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1231 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1232
1233 // pt trigger
1234 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1235 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1236 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1237 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1238 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1239 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1240 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1241 }
1242
1243 // pt associated
1244 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1245 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1246 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1247 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1248 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1249 }
1250
1251
1252 // ============================================================================================
1253 // the same for event mixing
1254 AliTHn *fHistPMix = bfMix->GetHistNp();
1255 AliTHn *fHistNMix = bfMix->GetHistNn();
1256 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1257 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1258 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1259 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1260
1261 // Psi_2
622473b9 1262 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1263 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1264 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1265 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1266 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1267 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1268
1269 // pt trigger
1270 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1271 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1272 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1273 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1274 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1275 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1276 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1277 }
1278
1279 // pt associated
1280 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1281 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1282 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1283 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1284 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1285 }
1286 // ============================================================================================
1287
1288
1289 //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)));
1290
1291 // Project into the wanted space (1st: analysis step, 2nd: axis)
1292 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1293 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1294 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1295 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1296 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1297 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1298
1299 // ============================================================================================
1300 // the same for event mixing
1301 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1302 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1303 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1304 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1305 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1306 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1307 // ============================================================================================
1308
1309 TH1D *gHistBalanceFunctionHistogram = 0x0;
1310 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1311
1312 hTemp1->Sumw2();
1313 hTemp2->Sumw2();
1314 hTemp3->Sumw2();
1315 hTemp4->Sumw2();
1316 hTemp1Mix->Sumw2();
1317 hTemp2Mix->Sumw2();
1318 hTemp3Mix->Sumw2();
1319 hTemp4Mix->Sumw2();
1320
1321 hTemp1->Scale(1./hTemp5->GetEntries());
1322 hTemp3->Scale(1./hTemp5->GetEntries());
1323 hTemp2->Scale(1./hTemp6->GetEntries());
1324 hTemp4->Scale(1./hTemp6->GetEntries());
1325
1326 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1327 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1328 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1329 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
1330
1331 hTemp1->Divide(hTemp1Mix);
1332 hTemp2->Divide(hTemp2Mix);
1333 hTemp3->Divide(hTemp3Mix);
1334 hTemp4->Divide(hTemp4Mix);
1335
1336 // now only project on one axis
1337 TH1D *h1DTemp1 = NULL;
1338 TH1D *h1DTemp2 = NULL;
1339 TH1D *h1DTemp3 = NULL;
1340 TH1D *h1DTemp4 = NULL;
1341 if(!bPhi){
1342 h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
1343 h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
1344 h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
1345 h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
1346 }
1347 else{
1348 h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
1349 h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
1350 h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
1351 h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
1352 }
1353
1354 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1355 gHistBalanceFunctionHistogram->Reset();
1356 if(!bPhi){
1357 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1358 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1359 }
1360 else{
1361 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1362 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1363 }
1364
1365 h1DTemp1->Add(h1DTemp3,-1.);
1366 h1DTemp2->Add(h1DTemp4,-1.);
1367
1368 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1369 gHistBalanceFunctionHistogram->Scale(0.5);
1370 }
1371
1372 return gHistBalanceFunctionHistogram;
1373}
1374
f0c5040c 1375
621329de 1376//____________________________________________________________________//
1377TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1378 Double_t psiMax,
1379 Double_t ptTriggerMin,
1380 Double_t ptTriggerMax,
1381 Double_t ptAssociatedMin,
47e040b5 1382 Double_t ptAssociatedMax) {
a38fd7f3 1383 //Returns the 2D correlation function for (+-) pairs
622473b9 1384
1385 // security checks
b0d53e32 1386 if(psiMin > psiMax-0.00001){
1387 AliError("psiMax <= psiMin");
1388 return NULL;
1389 }
1390 if(ptTriggerMin > ptTriggerMax-0.00001){
1391 AliError("ptTriggerMax <= ptTriggerMin");
1392 return NULL;
1393 }
1394 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1395 AliError("ptAssociatedMax <= ptAssociatedMin");
1396 return NULL;
1397 }
622473b9 1398
621329de 1399 // Psi_2: axis 0
622473b9 1400 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1401 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1402
6acdbcb2 1403 // pt trigger
1404 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1405 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1406 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1407 }
621329de 1408
6acdbcb2 1409 // pt associated
1410 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1411 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1412
1413 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1414 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1415
1416 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1417 //TCanvas *c1 = new TCanvas("c1","");
1418 //c1->cd();
1419 //if(!gHistTest){
1420 //AliError("Projection of fHistP = NULL");
1421 //return gHistTest;
1422 //}
1423 //else{
1424 //gHistTest->DrawCopy("colz");
1425 //}
1426
1427 //0:step, 1: Delta eta, 2: Delta phi
621329de 1428 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1429 if(!gHist){
1430 AliError("Projection of fHistPN = NULL");
1431 return gHist;
1432 }
1433
6acdbcb2 1434 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1435 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1436 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1437
6acdbcb2 1438 //TCanvas *c2 = new TCanvas("c2","");
1439 //c2->cd();
1440 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1441
47e040b5 1442 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1443 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1444
1445 //normalize to bin width
1446 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1447
1448 return gHist;
1449}
1450
1451//____________________________________________________________________//
621329de 1452TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1453 Double_t psiMax,
1454 Double_t ptTriggerMin,
1455 Double_t ptTriggerMax,
1456 Double_t ptAssociatedMin,
47e040b5 1457 Double_t ptAssociatedMax) {
a38fd7f3 1458 //Returns the 2D correlation function for (+-) pairs
622473b9 1459
1460 // security checks
b0d53e32 1461 if(psiMin > psiMax-0.00001){
1462 AliError("psiMax <= psiMin");
1463 return NULL;
1464 }
1465 if(ptTriggerMin > ptTriggerMax-0.00001){
1466 AliError("ptTriggerMax <= ptTriggerMin");
1467 return NULL;
1468 }
1469 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1470 AliError("ptAssociatedMax <= ptAssociatedMin");
1471 return NULL;
1472 }
622473b9 1473
621329de 1474 // Psi_2: axis 0
622473b9 1475 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1476 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
a38fd7f3 1477
6acdbcb2 1478 // pt trigger
1479 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1480 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1481 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1482 }
1483
1484 // pt associated
1485 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1486 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1487
1488 //0:step, 1: Delta eta, 2: Delta phi
621329de 1489 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1490 if(!gHist){
1491 AliError("Projection of fHistPN = NULL");
1492 return gHist;
1493 }
1494
a38fd7f3 1495 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1496 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 1497 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1498 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1499
1500 //normalize to bin width
1501 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1502
1503 return gHist;
1504}
1505
1506//____________________________________________________________________//
621329de 1507TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1508 Double_t psiMax,
1509 Double_t ptTriggerMin,
1510 Double_t ptTriggerMax,
1511 Double_t ptAssociatedMin,
47e040b5 1512 Double_t ptAssociatedMax) {
a38fd7f3 1513 //Returns the 2D correlation function for (+-) pairs
622473b9 1514
1515 // security checks
b0d53e32 1516 if(psiMin > psiMax-0.00001){
1517 AliError("psiMax <= psiMin");
1518 return NULL;
1519 }
1520 if(ptTriggerMin > ptTriggerMax-0.00001){
1521 AliError("ptTriggerMax <= ptTriggerMin");
1522 return NULL;
1523 }
1524 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1525 AliError("ptAssociatedMax <= ptAssociatedMin");
1526 return NULL;
1527 }
622473b9 1528
621329de 1529 // Psi_2: axis 0
622473b9 1530 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1531 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1532
1533 // pt trigger
1534 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1535 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1536 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1537 }
1538
1539 // pt associated
1540 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1541 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1542
6acdbcb2 1543 //0:step, 1: Delta eta, 2: Delta phi
621329de 1544 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1545 if(!gHist){
1546 AliError("Projection of fHistPN = NULL");
1547 return gHist;
1548 }
1549
a38fd7f3 1550 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1551 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 1552 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1553 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1554
1555 //normalize to bin width
1556 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1557
1558 return gHist;
1559}
1560
1561//____________________________________________________________________//
621329de 1562TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1563 Double_t psiMax,
1564 Double_t ptTriggerMin,
1565 Double_t ptTriggerMax,
1566 Double_t ptAssociatedMin,
47e040b5 1567 Double_t ptAssociatedMax) {
a38fd7f3 1568 //Returns the 2D correlation function for (+-) pairs
622473b9 1569
1570 // security checks
1571 if(psiMin > psiMax-0.00001){
1572 AliError("psiMax <= psiMin");
1573 return NULL;
1574 }
1575 if(ptTriggerMin > ptTriggerMax-0.00001){
1576 AliError("ptTriggerMax <= ptTriggerMin");
1577 return NULL;
1578 }
1579 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1580 AliError("ptAssociatedMax <= ptAssociatedMin");
1581 return NULL;
1582 }
1583
621329de 1584 // Psi_2: axis 0
622473b9 1585 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1586 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1587
1588 // pt trigger
1589 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1590 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1591 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1592 }
1593
1594 // pt associated
1595 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1596 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1597
6acdbcb2 1598 //0:step, 1: Delta eta, 2: Delta phi
621329de 1599 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 1600 if(!gHist){
1601 AliError("Projection of fHistPN = NULL");
1602 return gHist;
1603 }
1604
a38fd7f3 1605 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1606 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 1607 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1608 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1609
1610 //normalize to bin width
1611 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1612
1613 return gHist;
1614}
621329de 1615
7b610f27 1616//____________________________________________________________________//
1617TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1618 Double_t psiMax,
1619 Double_t ptTriggerMin,
1620 Double_t ptTriggerMax,
1621 Double_t ptAssociatedMin,
47e040b5 1622 Double_t ptAssociatedMax) {
7b610f27 1623 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 1624
1625 // security checks
b0d53e32 1626 if(psiMin > psiMax-0.00001){
1627 AliError("psiMax <= psiMin");
1628 return NULL;
1629 }
1630 if(ptTriggerMin > ptTriggerMax-0.00001){
1631 AliError("ptTriggerMax <= ptTriggerMin");
1632 return NULL;
1633 }
1634 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1635 AliError("ptAssociatedMax <= ptAssociatedMin");
1636 return NULL;
1637 }
622473b9 1638
7b610f27 1639 // Psi_2: axis 0
622473b9 1640 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1641 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1642 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1643 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1644 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1645 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
7b610f27 1646
1647 // pt trigger
1648 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1649 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1650 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1651 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1652 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1653 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1654 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 1655 }
1656
1657 // pt associated
1658 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1659 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1660 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1661 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1662 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 1663
1664 //0:step, 1: Delta eta, 2: Delta phi
1665 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1666 if(!gHistNN){
1667 AliError("Projection of fHistNN = NULL");
1668 return gHistNN;
1669 }
1670 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1671 if(!gHistPP){
1672 AliError("Projection of fHistPP = NULL");
1673 return gHistPP;
1674 }
1675 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1676 if(!gHistNP){
1677 AliError("Projection of fHistNP = NULL");
1678 return gHistNP;
1679 }
1680 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1681 if(!gHistPN){
1682 AliError("Projection of fHistPN = NULL");
1683 return gHistPN;
1684 }
1685
1686 // sum all 2 particle histograms
1687 gHistNN->Add(gHistPP);
1688 gHistNN->Add(gHistNP);
1689 gHistNN->Add(gHistPN);
1690
47e040b5 1691 // divide by sum of + and - triggers
1692 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 1693 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 1694
1695 //normalize to bin width
1696 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 1697
1698 return gHistNN;
1699}
1700
f2e8af26 1701//____________________________________________________________________//
1702
1703Bool_t AliBalancePsi::GetMomentsAnalytical(TH1D* gHist,
1704 Double_t &mean, Double_t &meanError,
1705 Double_t &sigma, Double_t &sigmaError,
1706 Double_t &skewness, Double_t &skewnessError,
1707 Double_t &kurtosis, Double_t &kurtosisError) {
1708 //
1709 // helper method to calculate the moments and errors of a TH1D anlytically
1710 //
1711
1712 Bool_t success = kFALSE;
1713 mean = -1.;
1714 meanError = -1.;
1715 sigma = -1.;
1716 sigmaError = -1.;
1717 skewness = -1.;
1718 skewnessError = -1.;
1719 kurtosis = -1.;
1720 kurtosisError = -1.;
1721
1722 if(gHist){
1723
1724 // ----------------------------------------------------------------------
1725 // basic parameters of histogram
1726
1727 Int_t fNumberOfBins = gHist->GetNbinsX();
1728 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
1729
1730
1731 // ----------------------------------------------------------------------
1732 // first calculate the mean
1733
1734 Double_t fWeightedAverage = 0.;
f0754617 1735 Double_t fNormalization = 0.;
f2e8af26 1736
1737 for(Int_t i = 1; i <= fNumberOfBins; i++) {
1738
1739 fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
1740 fNormalization += gHist->GetBinContent(i);
f2e8af26 1741 }
1742
1743 mean = fWeightedAverage / fNormalization;
1744
f2e8af26 1745 // ----------------------------------------------------------------------
1746 // then calculate the higher moments
1747
1748 Double_t fDelta = 0.;
1749 Double_t fDelta2 = 0.;
1750 Double_t fDelta3 = 0.;
1751 Double_t fDelta4 = 0.;
1752
1753 for(Int_t i = 1; i <= fNumberOfBins; i++) {
f0754617 1754 fDelta += gHist->GetBinContent(i) * (gHist->GetBinCenter(i) - mean);
1755 fDelta2 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
1756 fDelta3 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
1757 fDelta4 += gHist->GetBinContent(i) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
f2e8af26 1758 }
1759
f0754617 1760 sigma = TMath::Sqrt(fDelta2 / fNormalization);
f2e8af26 1761 skewness = fDelta3 / fNormalization / TMath::Power(sigma,3);
1762 kurtosis = fDelta4 / fNormalization / TMath::Power(sigma,4) - 3;
1763
f0754617 1764 // ----------------------------------------------------------------------
1765 // then calculate the higher moment errors?
1766 cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
1767 meanError = sigma / TMath::Sqrt(fNormalization); //sigma / TMath::Sqrt(gHist->GetEffectiveEntries()) would give the "ROOT" menaError
1768
f2e8af26 1769 success = kTRUE;
1770 }
1771
1772
1773 return success;
1774}
1775
1776
73609a48 1777//____________________________________________________________________//
1778Float_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) {
1779 //
1780 // calculates dphistar
1781 //
1782 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1783
1784 static const Double_t kPi = TMath::Pi();
1785
1786 // circularity
1787// if (dphistar > 2 * kPi)
1788// dphistar -= 2 * kPi;
1789// if (dphistar < -2 * kPi)
1790// dphistar += 2 * kPi;
1791
1792 if (dphistar > kPi)
1793 dphistar = kPi * 2 - dphistar;
1794 if (dphistar < -kPi)
1795 dphistar = -kPi * 2 - dphistar;
1796 if (dphistar > kPi) // might look funny but is needed
1797 dphistar = kPi * 2 - dphistar;
1798
1799 return dphistar;
1800}
1801
2922bbe3 1802
1803