]> 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
adjust output filename for LEGO
[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
b0d53e32 613 if(psiMin > psiMax-0.00001){
614 AliError("psiMax <= psiMin");
615 return NULL;
616 }
617 if(ptTriggerMin > ptTriggerMax-0.00001){
618 AliError("ptTriggerMax <= ptTriggerMin");
619 return NULL;
620 }
621 if(ptAssociatedMin > ptAssociatedMax-0.00001){
622 AliError("ptAssociatedMax <= ptAssociatedMin");
623 return NULL;
624 }
622473b9 625
9afe3098 626 // Psi_2
622473b9 627 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
628 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
629 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
630 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
631 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
632 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
9afe3098 633
6acdbcb2 634 // pt trigger
635 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 636 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
637 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
638 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
639 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
640 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
641 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 642 }
643
644 // pt associated
645 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 646 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
647 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
648 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
649 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 650 }
651
a38fd7f3 652 //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));
653
9afe3098 654 // Project into the wanted space (1st: analysis step, 2nd: axis)
655 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
656 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
657 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
658 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
659 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
660 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
661
662 TH1D *gHistBalanceFunctionHistogram = 0x0;
663 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
664 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
665 gHistBalanceFunctionHistogram->Reset();
666
667 switch(iVariablePair) {
621329de 668 case 1:
c8ad9635 669 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
670 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 671 break;
621329de 672 case 2:
c8ad9635 673 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
674 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 675 break;
9afe3098 676 default:
677 break;
678 }
0879e280 679
0879e280 680 hTemp1->Sumw2();
681 hTemp2->Sumw2();
682 hTemp3->Sumw2();
683 hTemp4->Sumw2();
a38fd7f3 684 hTemp1->Add(hTemp3,-1.);
0879e280 685 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 686 hTemp2->Add(hTemp4,-1.);
0879e280 687 hTemp2->Scale(1./hTemp6->GetEntries());
688 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 689 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 690
691 //normalize to bin width
692 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
0879e280 693 }
694
0879e280 695 return gHistBalanceFunctionHistogram;
696}
a38fd7f3 697
f0c5040c 698//____________________________________________________________________//
699TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
700 Int_t iVariablePair,
701 Double_t psiMin,
702 Double_t psiMax,
703 Double_t ptTriggerMin,
704 Double_t ptTriggerMax,
705 Double_t ptAssociatedMin,
706 Double_t ptAssociatedMax,
707 AliBalancePsi *bfMix) {
708 //Returns the BF histogram, extracted from the 6 AliTHn objects
709 //after dividing each correlation function by the Event Mixing one
710 //(private members) of the AliBalancePsi class.
711 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
712 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 713
714 // security checks
b0d53e32 715 if(psiMin > psiMax-0.00001){
716 AliError("psiMax <= psiMin");
717 return NULL;
718 }
719 if(ptTriggerMin > ptTriggerMax-0.00001){
720 AliError("ptTriggerMax <= ptTriggerMin");
721 return NULL;
722 }
723 if(ptAssociatedMin > ptAssociatedMax-0.00001){
724 AliError("ptAssociatedMax <= ptAssociatedMin");
725 return NULL;
726 }
622473b9 727
f0c5040c 728 // Psi_2
622473b9 729 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
730 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
731 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
732 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
733 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
734 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 735
736 // pt trigger
737 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 738 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
739 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
740 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
741 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
742 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
743 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 744 }
745
746 // pt associated
747 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 748 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
749 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
750 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
751 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 752 }
753
754 // ============================================================================================
755 // the same for event mixing
756 AliTHn *fHistPMix = bfMix->GetHistNp();
757 AliTHn *fHistNMix = bfMix->GetHistNn();
758 AliTHn *fHistPNMix = bfMix->GetHistNpn();
759 AliTHn *fHistNPMix = bfMix->GetHistNnp();
760 AliTHn *fHistPPMix = bfMix->GetHistNpp();
761 AliTHn *fHistNNMix = bfMix->GetHistNnn();
762
763 // Psi_2
622473b9 764 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
765 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
766 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
767 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
768 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
769 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 770
771 // pt trigger
772 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 773 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
774 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
775 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
776 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
777 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
778 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 779 }
780
781 // pt associated
782 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 783 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
784 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
785 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
786 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 787 }
788 // ============================================================================================
789
790 //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));
791
792 // Project into the wanted space (1st: analysis step, 2nd: axis)
793 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
794 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
795 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
796 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
797 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
798 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
799
800 // ============================================================================================
801 // the same for event mixing
f0fb4ac1 802 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
803 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
804 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
805 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
806 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
807 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
f0c5040c 808 // ============================================================================================
809
810 TH1D *gHistBalanceFunctionHistogram = 0x0;
811 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
812 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
813 gHistBalanceFunctionHistogram->Reset();
814
815 switch(iVariablePair) {
816 case 1:
817 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
818 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
819 break;
820 case 2:
821 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
822 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
823 break;
824 default:
825 break;
826 }
827
828 hTemp1->Sumw2();
829 hTemp2->Sumw2();
830 hTemp3->Sumw2();
831 hTemp4->Sumw2();
832 hTemp1Mix->Sumw2();
833 hTemp2Mix->Sumw2();
834 hTemp3Mix->Sumw2();
835 hTemp4Mix->Sumw2();
836
837 hTemp1->Scale(1./hTemp5->GetEntries());
838 hTemp3->Scale(1./hTemp5->GetEntries());
839 hTemp2->Scale(1./hTemp6->GetEntries());
840 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 841 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
842 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
843 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
844 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 845
846 hTemp1->Divide(hTemp1Mix);
847 hTemp2->Divide(hTemp2Mix);
848 hTemp3->Divide(hTemp3Mix);
849 hTemp4->Divide(hTemp4Mix);
850
851 hTemp1->Add(hTemp3,-1.);
852 hTemp2->Add(hTemp4,-1.);
853
854 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
855 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 856
857 //normalize to bin width
858 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 859 }
860
861 return gHistBalanceFunctionHistogram;
862}
863
a38fd7f3 864//____________________________________________________________________//
621329de 865TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 866 Double_t psiMax,
867 Double_t ptTriggerMin,
868 Double_t ptTriggerMax,
869 Double_t ptAssociatedMin,
870 Double_t ptAssociatedMax) {
621329de 871 //Returns the BF histogram in Delta eta vs Delta phi,
872 //extracted from the 6 AliTHn objects
873 //(private members) of the AliBalancePsi class.
874 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
875 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
622473b9 876
877 // security checks
b0d53e32 878 if(psiMin > psiMax-0.00001){
879 AliError("psiMax <= psiMin");
880 return NULL;
881 }
882 if(ptTriggerMin > ptTriggerMax-0.00001){
883 AliError("ptTriggerMax <= ptTriggerMin");
884 return NULL;
885 }
886 if(ptAssociatedMin > ptAssociatedMax-0.00001){
887 AliError("ptAssociatedMax <= ptAssociatedMin");
888 return NULL;
889 }
622473b9 890
621329de 891 TString histName = "gHistBalanceFunctionHistogram2D";
892
893 // Psi_2
622473b9 894 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
895 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
896 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
897 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
898 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
899 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 900
6acdbcb2 901 // pt trigger
902 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 903 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
904 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
905 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
906 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
907 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
908 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 909 }
910
911 // pt associated
912 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 913 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
914 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
915 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
916 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 917 }
918
919 //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 920
921 // Project into the wanted space (1st: analysis step, 2nd: axis)
922 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
923 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
924 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
925 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
926 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
927 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
928
929 TH2D *gHistBalanceFunctionHistogram = 0x0;
930 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
931 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
932 gHistBalanceFunctionHistogram->Reset();
c8ad9635 933 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
934 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
935 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 936
937 hTemp1->Sumw2();
938 hTemp2->Sumw2();
939 hTemp3->Sumw2();
940 hTemp4->Sumw2();
941 hTemp1->Add(hTemp3,-1.);
942 hTemp1->Scale(1./hTemp5->GetEntries());
943 hTemp2->Add(hTemp4,-1.);
944 hTemp2->Scale(1./hTemp6->GetEntries());
945 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
946 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 947
948 //normalize to bin width
949 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
621329de 950 }
951
952 return gHistBalanceFunctionHistogram;
953}
954
f0c5040c 955//____________________________________________________________________//
956TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
957 Double_t psiMax,
958 Double_t ptTriggerMin,
959 Double_t ptTriggerMax,
960 Double_t ptAssociatedMin,
961 Double_t ptAssociatedMax,
962 AliBalancePsi *bfMix) {
963 //Returns the BF histogram in Delta eta vs Delta phi,
964 //after dividing each correlation function by the Event Mixing one
965 //extracted from the 6 AliTHn objects
966 //(private members) of the AliBalancePsi class.
967 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
968 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
969 TString histName = "gHistBalanceFunctionHistogram2D";
970
971 if(!bfMix){
972 AliError("balance function object for event mixing not available");
973 return NULL;
974 }
975
622473b9 976 // security checks
b0d53e32 977 if(psiMin > psiMax-0.00001){
978 AliError("psiMax <= psiMin");
979 return NULL;
980 }
981 if(ptTriggerMin > ptTriggerMax-0.00001){
982 AliError("ptTriggerMax <= ptTriggerMin");
983 return NULL;
984 }
985 if(ptAssociatedMin > ptAssociatedMax-0.00001){
986 AliError("ptAssociatedMax <= ptAssociatedMin");
987 return NULL;
988 }
622473b9 989
f0c5040c 990 // Psi_2
622473b9 991 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
992 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
993 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
994 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
995 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
996 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 997
998 // pt trigger
999 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1000 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1001 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1002 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1003 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1004 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1005 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1006 }
1007
1008 // pt associated
1009 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1010 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1011 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1012 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1013 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1014 }
1015
1016
1017 // ============================================================================================
1018 // the same for event mixing
1019 AliTHn *fHistPMix = bfMix->GetHistNp();
1020 AliTHn *fHistNMix = bfMix->GetHistNn();
1021 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1022 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1023 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1024 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1025
1026 // Psi_2
622473b9 1027 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1028 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1029 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1030 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1031 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1032 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f0c5040c 1033
1034 // pt trigger
1035 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1036 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1037 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1038 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1039 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1040 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1041 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f0c5040c 1042 }
1043
1044 // pt associated
1045 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1046 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1047 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1048 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1049 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f0c5040c 1050 }
1051 // ============================================================================================
1052
1053
1054 //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)));
1055
1056 // Project into the wanted space (1st: analysis step, 2nd: axis)
1057 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1058 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1059 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1060 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1061 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1062 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1063
1064 // ============================================================================================
1065 // the same for event mixing
1066 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1067 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1068 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1069 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
f0fb4ac1 1070 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1071 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
f0c5040c 1072 // ============================================================================================
1073
1074 TH2D *gHistBalanceFunctionHistogram = 0x0;
1075 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1076 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1077 gHistBalanceFunctionHistogram->Reset();
1078 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1079 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1080 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1081
1082 hTemp1->Sumw2();
1083 hTemp2->Sumw2();
1084 hTemp3->Sumw2();
1085 hTemp4->Sumw2();
1086 hTemp1Mix->Sumw2();
1087 hTemp2Mix->Sumw2();
1088 hTemp3Mix->Sumw2();
1089 hTemp4Mix->Sumw2();
1090
1091 hTemp1->Scale(1./hTemp5->GetEntries());
1092 hTemp3->Scale(1./hTemp5->GetEntries());
1093 hTemp2->Scale(1./hTemp6->GetEntries());
1094 hTemp4->Scale(1./hTemp6->GetEntries());
f0fb4ac1 1095 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1096 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1097 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1098 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
f0c5040c 1099
1100 hTemp1->Divide(hTemp1Mix);
1101 hTemp2->Divide(hTemp2Mix);
1102 hTemp3->Divide(hTemp3Mix);
1103 hTemp4->Divide(hTemp4Mix);
1104
1105 hTemp1->Add(hTemp3,-1.);
1106 hTemp2->Add(hTemp4,-1.);
1107
1108 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1109 gHistBalanceFunctionHistogram->Scale(0.5);
faa03677 1110
1111 //normalize to bin width
1112 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
f0c5040c 1113 }
1114
1115 return gHistBalanceFunctionHistogram;
1116}
1117
f2e8af26 1118//____________________________________________________________________//
1119TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1120 Double_t psiMin,
1121 Double_t psiMax,
1122 Double_t ptTriggerMin,
1123 Double_t ptTriggerMax,
1124 Double_t ptAssociatedMin,
1125 Double_t ptAssociatedMax,
1126 AliBalancePsi *bfMix) {
1127 //Returns the BF histogram in Delta eta OR Delta phi,
1128 //after dividing each correlation function by the Event Mixing one
1129 // (But the division is done here in 2D, this was basically done to check the Integral)
1130 //extracted from the 6 AliTHn objects
1131 //(private members) of the AliBalancePsi class.
1132 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1133 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1134 TString histName = "gHistBalanceFunctionHistogram1D";
1135
1136 if(!bfMix){
1137 AliError("balance function object for event mixing not available");
1138 return NULL;
1139 }
1140
622473b9 1141 // security checks
b0d53e32 1142 if(psiMin > psiMax-0.00001){
1143 AliError("psiMax <= psiMin");
1144 return NULL;
1145 }
1146 if(ptTriggerMin > ptTriggerMax-0.00001){
1147 AliError("ptTriggerMax <= ptTriggerMin");
1148 return NULL;
1149 }
1150 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1151 AliError("ptAssociatedMax <= ptAssociatedMin");
1152 return NULL;
1153 }
622473b9 1154
f2e8af26 1155 // Psi_2
622473b9 1156 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1157 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1158 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1159 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1160 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1161 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1162
1163 // pt trigger
1164 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1165 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1166 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1167 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1168 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1169 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1170 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1171 }
1172
1173 // pt associated
1174 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1175 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1176 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1177 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1178 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1179 }
1180
1181
1182 // ============================================================================================
1183 // the same for event mixing
1184 AliTHn *fHistPMix = bfMix->GetHistNp();
1185 AliTHn *fHistNMix = bfMix->GetHistNn();
1186 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1187 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1188 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1189 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1190
1191 // Psi_2
622473b9 1192 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1193 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1194 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1195 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1196 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1197 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
f2e8af26 1198
1199 // pt trigger
1200 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1201 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1202 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1203 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1204 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1205 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1206 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
f2e8af26 1207 }
1208
1209 // pt associated
1210 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
622473b9 1211 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1212 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1213 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1214 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
f2e8af26 1215 }
1216 // ============================================================================================
1217
1218
1219 //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)));
1220
1221 // Project into the wanted space (1st: analysis step, 2nd: axis)
1222 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1223 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1224 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1225 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1226 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1227 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1228
1229 // ============================================================================================
1230 // the same for event mixing
1231 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1232 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1233 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1234 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1235 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1236 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1237 // ============================================================================================
1238
1239 TH1D *gHistBalanceFunctionHistogram = 0x0;
1240 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1241
1242 hTemp1->Sumw2();
1243 hTemp2->Sumw2();
1244 hTemp3->Sumw2();
1245 hTemp4->Sumw2();
1246 hTemp1Mix->Sumw2();
1247 hTemp2Mix->Sumw2();
1248 hTemp3Mix->Sumw2();
1249 hTemp4Mix->Sumw2();
1250
1251 hTemp1->Scale(1./hTemp5->GetEntries());
1252 hTemp3->Scale(1./hTemp5->GetEntries());
1253 hTemp2->Scale(1./hTemp6->GetEntries());
1254 hTemp4->Scale(1./hTemp6->GetEntries());
1255
1256 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
1257 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
1258 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
1259 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
1260
1261 hTemp1->Divide(hTemp1Mix);
1262 hTemp2->Divide(hTemp2Mix);
1263 hTemp3->Divide(hTemp3Mix);
1264 hTemp4->Divide(hTemp4Mix);
1265
1266 // now only project on one axis
1267 TH1D *h1DTemp1 = NULL;
1268 TH1D *h1DTemp2 = NULL;
1269 TH1D *h1DTemp3 = NULL;
1270 TH1D *h1DTemp4 = NULL;
1271 if(!bPhi){
1272 h1DTemp1 = (TH1D*)hTemp1->ProjectionX(Form("%s_projX",hTemp1->GetName()));
1273 h1DTemp2 = (TH1D*)hTemp2->ProjectionX(Form("%s_projX",hTemp2->GetName()));
1274 h1DTemp3 = (TH1D*)hTemp3->ProjectionX(Form("%s_projX",hTemp3->GetName()));
1275 h1DTemp4 = (TH1D*)hTemp4->ProjectionX(Form("%s_projX",hTemp4->GetName()));
1276 }
1277 else{
1278 h1DTemp1 = (TH1D*)hTemp1->ProjectionY(Form("%s_projX",hTemp1->GetName()));
1279 h1DTemp2 = (TH1D*)hTemp2->ProjectionY(Form("%s_projX",hTemp2->GetName()));
1280 h1DTemp3 = (TH1D*)hTemp3->ProjectionY(Form("%s_projX",hTemp3->GetName()));
1281 h1DTemp4 = (TH1D*)hTemp4->ProjectionY(Form("%s_projX",hTemp4->GetName()));
1282 }
1283
1284 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1285 gHistBalanceFunctionHistogram->Reset();
1286 if(!bPhi){
1287 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1288 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1289 }
1290 else{
1291 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1292 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1293 }
1294
1295 h1DTemp1->Add(h1DTemp3,-1.);
1296 h1DTemp2->Add(h1DTemp4,-1.);
1297
1298 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1299 gHistBalanceFunctionHistogram->Scale(0.5);
1300 }
1301
1302 return gHistBalanceFunctionHistogram;
1303}
1304
f0c5040c 1305
621329de 1306//____________________________________________________________________//
1307TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 1308 Double_t psiMax,
1309 Double_t ptTriggerMin,
1310 Double_t ptTriggerMax,
1311 Double_t ptAssociatedMin,
47e040b5 1312 Double_t ptAssociatedMax) {
a38fd7f3 1313 //Returns the 2D correlation function for (+-) pairs
622473b9 1314
1315 // security checks
b0d53e32 1316 if(psiMin > psiMax-0.00001){
1317 AliError("psiMax <= psiMin");
1318 return NULL;
1319 }
1320 if(ptTriggerMin > ptTriggerMax-0.00001){
1321 AliError("ptTriggerMax <= ptTriggerMin");
1322 return NULL;
1323 }
1324 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1325 AliError("ptAssociatedMax <= ptAssociatedMin");
1326 return NULL;
1327 }
622473b9 1328
621329de 1329 // Psi_2: axis 0
622473b9 1330 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1331 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
621329de 1332
6acdbcb2 1333 // pt trigger
1334 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1335 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1336 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
5ef8eaa7 1337 }
621329de 1338
6acdbcb2 1339 // pt associated
1340 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1341 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1342
1343 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1344 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
1345
1346 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
1347 //TCanvas *c1 = new TCanvas("c1","");
1348 //c1->cd();
1349 //if(!gHistTest){
1350 //AliError("Projection of fHistP = NULL");
1351 //return gHistTest;
1352 //}
1353 //else{
1354 //gHistTest->DrawCopy("colz");
1355 //}
1356
1357 //0:step, 1: Delta eta, 2: Delta phi
621329de 1358 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 1359 if(!gHist){
1360 AliError("Projection of fHistPN = NULL");
1361 return gHist;
1362 }
1363
6acdbcb2 1364 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
1365 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
1366 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 1367
6acdbcb2 1368 //TCanvas *c2 = new TCanvas("c2","");
1369 //c2->cd();
1370 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 1371
47e040b5 1372 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1373 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1374
1375 //normalize to bin width
1376 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1377
1378 return gHist;
1379}
1380
1381//____________________________________________________________________//
621329de 1382TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 1383 Double_t psiMax,
1384 Double_t ptTriggerMin,
1385 Double_t ptTriggerMax,
1386 Double_t ptAssociatedMin,
47e040b5 1387 Double_t ptAssociatedMax) {
a38fd7f3 1388 //Returns the 2D correlation function for (+-) pairs
622473b9 1389
1390 // security checks
b0d53e32 1391 if(psiMin > psiMax-0.00001){
1392 AliError("psiMax <= psiMin");
1393 return NULL;
1394 }
1395 if(ptTriggerMin > ptTriggerMax-0.00001){
1396 AliError("ptTriggerMax <= ptTriggerMin");
1397 return NULL;
1398 }
1399 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1400 AliError("ptAssociatedMax <= ptAssociatedMin");
1401 return NULL;
1402 }
622473b9 1403
621329de 1404 // Psi_2: axis 0
622473b9 1405 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1406 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
a38fd7f3 1407
6acdbcb2 1408 // pt trigger
1409 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1410 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1411 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1412 }
1413
1414 // pt associated
1415 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1416 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
6acdbcb2 1417
1418 //0:step, 1: Delta eta, 2: Delta phi
621329de 1419 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 1420 if(!gHist){
1421 AliError("Projection of fHistPN = NULL");
1422 return gHist;
1423 }
1424
a38fd7f3 1425 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1426 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
47e040b5 1427 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1428 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1429
1430 //normalize to bin width
1431 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1432
1433 return gHist;
1434}
1435
1436//____________________________________________________________________//
621329de 1437TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1438 Double_t psiMax,
1439 Double_t ptTriggerMin,
1440 Double_t ptTriggerMax,
1441 Double_t ptAssociatedMin,
47e040b5 1442 Double_t ptAssociatedMax) {
a38fd7f3 1443 //Returns the 2D correlation function for (+-) pairs
622473b9 1444
1445 // security checks
b0d53e32 1446 if(psiMin > psiMax-0.00001){
1447 AliError("psiMax <= psiMin");
1448 return NULL;
1449 }
1450 if(ptTriggerMin > ptTriggerMax-0.00001){
1451 AliError("ptTriggerMax <= ptTriggerMin");
1452 return NULL;
1453 }
1454 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1455 AliError("ptAssociatedMax <= ptAssociatedMin");
1456 return NULL;
1457 }
622473b9 1458
621329de 1459 // Psi_2: axis 0
622473b9 1460 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1461 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1462
1463 // pt trigger
1464 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1465 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1466 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1467 }
1468
1469 // pt associated
1470 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1471 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1472
6acdbcb2 1473 //0:step, 1: Delta eta, 2: Delta phi
621329de 1474 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1475 if(!gHist){
1476 AliError("Projection of fHistPN = NULL");
1477 return gHist;
1478 }
1479
a38fd7f3 1480 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1481 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
47e040b5 1482 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
621329de 1483 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
faa03677 1484
1485 //normalize to bin width
1486 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1487
1488 return gHist;
1489}
1490
1491//____________________________________________________________________//
621329de 1492TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1493 Double_t psiMax,
1494 Double_t ptTriggerMin,
1495 Double_t ptTriggerMax,
1496 Double_t ptAssociatedMin,
47e040b5 1497 Double_t ptAssociatedMax) {
a38fd7f3 1498 //Returns the 2D correlation function for (+-) pairs
622473b9 1499
1500 // security checks
1501 if(psiMin > psiMax-0.00001){
1502 AliError("psiMax <= psiMin");
1503 return NULL;
1504 }
1505 if(ptTriggerMin > ptTriggerMax-0.00001){
1506 AliError("ptTriggerMax <= ptTriggerMin");
1507 return NULL;
1508 }
1509 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1510 AliError("ptAssociatedMax <= ptAssociatedMin");
1511 return NULL;
1512 }
1513
621329de 1514 // Psi_2: axis 0
622473b9 1515 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1516 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
6acdbcb2 1517
1518 // pt trigger
1519 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1520 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1521 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
6acdbcb2 1522 }
1523
1524 // pt associated
1525 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1526 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
a38fd7f3 1527
6acdbcb2 1528 //0:step, 1: Delta eta, 2: Delta phi
621329de 1529 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 1530 if(!gHist){
1531 AliError("Projection of fHistPN = NULL");
1532 return gHist;
1533 }
1534
a38fd7f3 1535 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1536 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
47e040b5 1537 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
621329de 1538 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
faa03677 1539
1540 //normalize to bin width
1541 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
a38fd7f3 1542
1543 return gHist;
1544}
621329de 1545
7b610f27 1546//____________________________________________________________________//
1547TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
1548 Double_t psiMax,
1549 Double_t ptTriggerMin,
1550 Double_t ptTriggerMax,
1551 Double_t ptAssociatedMin,
47e040b5 1552 Double_t ptAssociatedMax) {
7b610f27 1553 //Returns the 2D correlation function for the sum of all charge combination pairs
622473b9 1554
1555 // security checks
b0d53e32 1556 if(psiMin > psiMax-0.00001){
1557 AliError("psiMax <= psiMin");
1558 return NULL;
1559 }
1560 if(ptTriggerMin > ptTriggerMax-0.00001){
1561 AliError("ptTriggerMax <= ptTriggerMin");
1562 return NULL;
1563 }
1564 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1565 AliError("ptAssociatedMax <= ptAssociatedMin");
1566 return NULL;
1567 }
622473b9 1568
7b610f27 1569 // Psi_2: axis 0
622473b9 1570 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1571 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1572 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1573 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1574 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1575 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
7b610f27 1576
1577 // pt trigger
1578 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
622473b9 1579 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1580 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1581 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1582 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1583 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1584 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
7b610f27 1585 }
1586
1587 // pt associated
1588 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
622473b9 1589 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1590 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1591 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1592 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
7b610f27 1593
1594 //0:step, 1: Delta eta, 2: Delta phi
1595 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
1596 if(!gHistNN){
1597 AliError("Projection of fHistNN = NULL");
1598 return gHistNN;
1599 }
1600 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
1601 if(!gHistPP){
1602 AliError("Projection of fHistPP = NULL");
1603 return gHistPP;
1604 }
1605 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
1606 if(!gHistNP){
1607 AliError("Projection of fHistNP = NULL");
1608 return gHistNP;
1609 }
1610 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
1611 if(!gHistPN){
1612 AliError("Projection of fHistPN = NULL");
1613 return gHistPN;
1614 }
1615
1616 // sum all 2 particle histograms
1617 gHistNN->Add(gHistPP);
1618 gHistNN->Add(gHistNP);
1619 gHistNN->Add(gHistPN);
1620
47e040b5 1621 // divide by sum of + and - triggers
1622 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0 && (Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
7b610f27 1623 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries() + fHistN->Project(0,1)->GetEntries()));
faa03677 1624
1625 //normalize to bin width
1626 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
7b610f27 1627
1628 return gHistNN;
1629}
1630
f2e8af26 1631//____________________________________________________________________//
1632
1633Bool_t AliBalancePsi::GetMomentsAnalytical(TH1D* gHist,
1634 Double_t &mean, Double_t &meanError,
1635 Double_t &sigma, Double_t &sigmaError,
1636 Double_t &skewness, Double_t &skewnessError,
1637 Double_t &kurtosis, Double_t &kurtosisError) {
1638 //
1639 // helper method to calculate the moments and errors of a TH1D anlytically
1640 //
1641
1642 Bool_t success = kFALSE;
1643 mean = -1.;
1644 meanError = -1.;
1645 sigma = -1.;
1646 sigmaError = -1.;
1647 skewness = -1.;
1648 skewnessError = -1.;
1649 kurtosis = -1.;
1650 kurtosisError = -1.;
1651
1652 if(gHist){
1653
1654 // ----------------------------------------------------------------------
1655 // basic parameters of histogram
1656
1657 Int_t fNumberOfBins = gHist->GetNbinsX();
1658 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
1659
1660
1661 // ----------------------------------------------------------------------
1662 // first calculate the mean
1663
1664 Double_t fWeightedAverage = 0.;
1665 Double_t fNormalization = 0.;
1666
1667 for(Int_t i = 1; i <= fNumberOfBins; i++) {
1668
1669 fWeightedAverage += gHist->GetBinContent(i) * gHist->GetBinCenter(i);
1670 fNormalization += gHist->GetBinContent(i);
1671
1672 }
1673
1674 mean = fWeightedAverage / fNormalization;
1675
1676
1677 // ----------------------------------------------------------------------
1678 // then calculate the higher moments
1679
1680 Double_t fDelta = 0.;
1681 Double_t fDelta2 = 0.;
1682 Double_t fDelta3 = 0.;
1683 Double_t fDelta4 = 0.;
1684
1685 for(Int_t i = 1; i <= fNumberOfBins; i++) {
1686 fDelta += gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean;
1687 fDelta2 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),2);
1688 fDelta3 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),3);
1689 fDelta4 += TMath::Power((gHist->GetBinContent(i) * gHist->GetBinCenter(i) - mean),4);
1690 }
1691
1692 sigma = fDelta2 / fNormalization;
1693 skewness = fDelta3 / fNormalization / TMath::Power(sigma,3);
1694 kurtosis = fDelta4 / fNormalization / TMath::Power(sigma,4) - 3;
1695
1696 success = kTRUE;
1697 }
1698
1699
1700 return success;
1701}
1702
1703
73609a48 1704//____________________________________________________________________//
1705Float_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) {
1706 //
1707 // calculates dphistar
1708 //
1709 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1710
1711 static const Double_t kPi = TMath::Pi();
1712
1713 // circularity
1714// if (dphistar > 2 * kPi)
1715// dphistar -= 2 * kPi;
1716// if (dphistar < -2 * kPi)
1717// dphistar += 2 * kPi;
1718
1719 if (dphistar > kPi)
1720 dphistar = kPi * 2 - dphistar;
1721 if (dphistar < -kPi)
1722 dphistar = -kPi * 2 - dphistar;
1723 if (dphistar > kPi) // might look funny but is needed
1724 dphistar = kPi * 2 - dphistar;
1725
1726 return dphistar;
1727}
1728