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