]> 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
added ZNA histo (from Chiara)
[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"
0879e280 40
41#include "AliBalancePsi.h"
42
43ClassImp(AliBalancePsi)
44
45//____________________________________________________________________//
46AliBalancePsi::AliBalancePsi() :
47 TObject(),
9fd4b54e 48 fShuffle(kFALSE),
0879e280 49 fAnalysisLevel("ESD"),
50 fAnalyzedEvents(0) ,
51 fCentralityId(0) ,
52 fCentStart(0.),
53 fCentStop(0.),
9afe3098 54 fHistP(0),
55 fHistN(0),
56 fHistPN(0),
57 fHistNP(0),
58 fHistPP(0),
59 fHistNN(0),
73609a48 60 fHistHBTbefore(0),
61 fHistHBTafter(0),
62 fHistConversionbefore(0),
63 fHistConversionafter(0),
c8ad9635 64 fHistPsiMinusPhi(0),
73609a48 65 fPsiInterval(15.),
66 fHBTCut(kFALSE),
67 fConversionCut(kFALSE) {
0879e280 68 // Default constructor
0879e280 69}
70
0879e280 71//____________________________________________________________________//
72AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
9fd4b54e 73 TObject(balance), fShuffle(balance.fShuffle),
0879e280 74 fAnalysisLevel(balance.fAnalysisLevel),
75 fAnalyzedEvents(balance.fAnalyzedEvents),
76 fCentralityId(balance.fCentralityId),
77 fCentStart(balance.fCentStart),
78 fCentStop(balance.fCentStop),
9afe3098 79 fHistP(balance.fHistP),
80 fHistN(balance.fHistN),
81 fHistPN(balance.fHistPN),
82 fHistNP(balance.fHistNP),
83 fHistPP(balance.fHistPP),
84 fHistNN(balance.fHistNN),
73609a48 85 fHistHBTbefore(balance.fHistHBTbefore),
86 fHistHBTafter(balance.fHistHBTafter),
87 fHistConversionbefore(balance.fHistConversionbefore),
88 fHistConversionafter(balance.fHistConversionafter),
c8ad9635 89 fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
73609a48 90 fPsiInterval(balance.fPsiInterval),
91 fHBTCut(balance.fHBTCut),
92 fConversionCut(balance.fConversionCut) {
0879e280 93 //copy constructor
0879e280 94}
95
96//____________________________________________________________________//
97AliBalancePsi::~AliBalancePsi() {
98 // Destructor
9afe3098 99 delete fHistP;
100 delete fHistN;
101 delete fHistPN;
102 delete fHistNP;
103 delete fHistPP;
104 delete fHistNN;
73609a48 105
106 delete fHistHBTbefore;
107 delete fHistHBTafter;
108 delete fHistConversionbefore;
109 delete fHistConversionafter;
c8ad9635 110 delete fHistPsiMinusPhi;
0879e280 111}
112
113//____________________________________________________________________//
9afe3098 114void AliBalancePsi::InitHistograms() {
115 // single particle histograms
116 Int_t anaSteps = 1; // analysis steps
9fd4b54e 117 Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
118 Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
119 TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
9afe3098 120
121 // two particle histograms
9fd4b54e 122 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
123 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
124 TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
9afe3098 125
126 //centrality
621329de 127 /*const Int_t kNCentralityBins = 9;
9afe3098 128 Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
9afe3098 129 iBinSingle[0] = kNCentralityBins;
130 dBinsSingle[0] = centralityBins;
131 axisTitleSingle[0] = "Centrality percentile [%]";
132 iBinPair[0] = kNCentralityBins;
133 dBinsPair[0] = centralityBins;
621329de 134 axisTitlePair[0] = "Centrality percentile [%]"; */
9afe3098 135
c8ad9635 136 //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)
6acdbcb2 137 const Int_t kNPsi2Bins = 4;
138 Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
621329de 139 iBinSingle[0] = kNPsi2Bins;
140 dBinsSingle[0] = psi2Bins;
c8ad9635 141 axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
621329de 142 iBinPair[0] = kNPsi2Bins;
143 dBinsPair[0] = psi2Bins;
c8ad9635 144 axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
9afe3098 145
621329de 146 // delta eta
147 const Int_t kNDeltaEtaBins = 80;
9afe3098 148 Double_t deltaEtaBins[kNDeltaEtaBins+1];
149 for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
621329de 150 deltaEtaBins[i] = -2.0 + i * 0.05;
151 iBinPair[1] = kNDeltaEtaBins;
152 dBinsPair[1] = deltaEtaBins;
c8ad9635 153 axisTitlePair[1] = "#Delta#eta";
621329de 154
155 // delta phi
9afe3098 156 const Int_t kNDeltaPhiBins = 72;
157 Double_t deltaPhiBins[kNDeltaPhiBins+1];
158 for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
c8ad9635 159 //deltaPhiBins[i] = -180.0 + i * 5.;
160 deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.;
9afe3098 161 }
621329de 162 iBinPair[2] = kNDeltaPhiBins;
163 dBinsPair[2] = deltaPhiBins;
c8ad9635 164 axisTitlePair[2] = "#Delta#varphi (rad)";
9afe3098 165
a38fd7f3 166 // pt(trigger-associated)
8795e993 167 const Int_t kNPtBins = 16;
168 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.};
169 //for(Int_t i = 0; i < kNPtBins+1; i++){
170 //ptBins[i] = 0.2 + i * 0.5;
171 //}
621329de 172 iBinSingle[1] = kNPtBins;
173 dBinsSingle[1] = ptBins;
c8ad9635 174 axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
621329de 175
176 iBinPair[3] = kNPtBins;
177 dBinsPair[3] = ptBins;
c8ad9635 178 axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
a38fd7f3 179
180 iBinPair[4] = kNPtBins;
181 dBinsPair[4] = ptBins;
c8ad9635 182 axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
9afe3098 183
184 TString histName;
185 //+ triggered particles
186 histName = "fHistP";
9fd4b54e 187 if(fShuffle) histName.Append("_shuffle");
9afe3098 188 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 189 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
190 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 191 fHistP->SetBinLimits(j, dBinsSingle[j]);
192 fHistP->SetVarTitle(j, axisTitleSingle[j]);
0879e280 193 }
9afe3098 194
195 //- triggered particles
196 histName = "fHistN";
9fd4b54e 197 if(fShuffle) histName.Append("_shuffle");
9afe3098 198 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 199 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
200 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 201 fHistN->SetBinLimits(j, dBinsSingle[j]);
202 fHistN->SetVarTitle(j, axisTitleSingle[j]);
0879e280 203 }
9afe3098 204
205 //+- pairs
206 histName = "fHistPN";
9fd4b54e 207 if(fShuffle) histName.Append("_shuffle");
9afe3098 208 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 209 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
210 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 211 fHistPN->SetBinLimits(j, dBinsPair[j]);
212 fHistPN->SetVarTitle(j, axisTitlePair[j]);
0879e280 213 }
0879e280 214
9afe3098 215 //-+ pairs
216 histName = "fHistNP";
9fd4b54e 217 if(fShuffle) histName.Append("_shuffle");
9afe3098 218 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 219 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
220 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 221 fHistNP->SetBinLimits(j, dBinsPair[j]);
222 fHistNP->SetVarTitle(j, axisTitlePair[j]);
0879e280 223 }
0879e280 224
9afe3098 225 //++ pairs
226 histName = "fHistPP";
9fd4b54e 227 if(fShuffle) histName.Append("_shuffle");
9afe3098 228 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 229 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
230 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 231 fHistPP->SetBinLimits(j, dBinsPair[j]);
232 fHistPP->SetVarTitle(j, axisTitlePair[j]);
233 }
234
235 //-- pairs
236 histName = "fHistNN";
9fd4b54e 237 if(fShuffle) histName.Append("_shuffle");
9afe3098 238 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 239 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
240 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 241 fHistNN->SetBinLimits(j, dBinsPair[j]);
242 fHistNN->SetVarTitle(j, axisTitlePair[j]);
0879e280 243 }
f06d59b3 244 AliInfo("Finished setting up the AliTHn");
73609a48 245
246 // QA histograms
c8ad9635 247 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
248 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
249 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
250 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
251 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
0879e280 252}
253
254//____________________________________________________________________//
f06d59b3 255void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
256 TObjArray *particles,
73609a48 257 TObjArray *particlesMixed,
258 Float_t bSign) {
0879e280 259 // Calculates the balance function
260 fAnalyzedEvents++;
0879e280 261
262 // Initialize histograms if not done yet
9afe3098 263 if(!fHistPN){
0879e280 264 AliWarning("Histograms not yet initialized! --> Will be done now");
265 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
266 InitHistograms();
267 }
268
9fd4b54e 269 Double_t trackVariablesSingle[kTrackVariablesSingle];
270 Double_t trackVariablesPair[kTrackVariablesPair];
f06d59b3 271
272 if (!particles){
273 AliWarning("particles TObjArray is NULL pointer --> return");
274 return;
275 }
9afe3098 276
f06d59b3 277 // define end of particle loops
278 Int_t iMax = particles->GetEntriesFast();
279 Int_t jMax = iMax;
280 if (particlesMixed)
281 jMax = particlesMixed->GetEntriesFast();
282
283 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
284 TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
285
286 TArrayF secondEta(jMax);
287 TArrayF secondPhi(jMax);
288 TArrayF secondPt(jMax);
289
290 for (Int_t i=0; i<jMax; i++){
291 secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
292 secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
293 secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
294 }
295
296 // 1st particle loop
73609a48 297 for (Int_t i = 0; i < iMax; i++) {
6acdbcb2 298 AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
299
300 // some optimization
301 Float_t firstEta = firstParticle->Eta();
302 Float_t firstPhi = firstParticle->Phi();
303 Float_t firstPt = firstParticle->Pt();
304
305 // Event plane (determine psi bin)
306 Double_t gPsiMinusPhi = 0.;
307 Double_t gPsiMinusPhiBin = -10.;
308 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
309 //in-plane
c8ad9635 310 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
311 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
6acdbcb2 312 gPsiMinusPhiBin = 0.0;
313 //intermediate
c8ad9635 314 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
315 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
316 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
317 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
6acdbcb2 318 gPsiMinusPhiBin = 1.0;
319 //out of plane
c8ad9635 320 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
321 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
6acdbcb2 322 gPsiMinusPhiBin = 2.0;
323 //everything else
324 else
325 gPsiMinusPhiBin = 3.0;
326
c8ad9635 327 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
73609a48 328
329 Short_t charge1 = (Short_t) firstParticle->Charge();
6acdbcb2 330
331 trackVariablesSingle[0] = gPsiMinusPhiBin;
332 trackVariablesSingle[1] = firstPt;
333
334 //fill single particle histograms
73609a48 335 if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,1.);
336 else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,1.);
6acdbcb2 337
338 // 2nd particle loop (only for j < i for non double counting in the same pT region)
339 // --> SAME pT region for trigger and assoc: NO double counting with this
340 // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
341 // --> can be handled afterwards by using assoc. as trigger as well ?!
a86235ae 342 for(Int_t j = 0; j < iMax; j++) {
a86235ae 343 if (particlesMixed && j == jMax-1 ) // if the mixed track number is smaller than the main event one
6acdbcb2 344 break;
a86235ae 345
346 if(j == i) continue; // no auto correlations
9afe3098 347
6acdbcb2 348 AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
6acdbcb2 349 Short_t charge2 = (Short_t) secondParticle->Charge();
9afe3098 350
6acdbcb2 351 trackVariablesPair[0] = gPsiMinusPhiBin;
352 trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
353 trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
c8ad9635 354 //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
355 //trackVariablesPair[2] -= 360.;
356 //if (trackVariablesPair[2] < - 180.)
357 //trackVariablesPair[2] += 360.;
358 if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
359 trackVariablesPair[2] -= 2.*TMath::Pi();
360 if (trackVariablesPair[2] < - TMath::Pi())
361 trackVariablesPair[2] += 2.*TMath::Pi();
362 if (trackVariablesPair[2] < - TMath::Pi()/2.)
363 trackVariablesPair[2] += 2.*TMath::Pi();
9afe3098 364
6acdbcb2 365 trackVariablesPair[3] = firstPt; // pt trigger
366 trackVariablesPair[4] = secondPt[j]; // pt
367 // trackVariablesPair[5] = fCentrality; // centrality
73609a48 368
369 // HBT like cut
370 if(fHBTCut && charge1 * charge2 > 0){
371 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
372 // continue;
373
374 Double_t deta = firstEta - secondEta[j];
375 Double_t dphi = firstPhi - secondPhi[j];
376 // VERSION 2 (Taken from DPhiCorrelations)
377 // the variables & cuthave been developed by the HBT group
378 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
379 fHistHBTbefore->Fill(deta,dphi);
380
381 // optimization
382 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
383 {
384 // phi in rad
c8ad9635 385 //Float_t phi1rad = firstPhi*TMath::DegToRad();
386 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
387 Float_t phi1rad = firstPhi;
388 Float_t phi2rad = secondPhi[j];
73609a48 389
390 // check first boundaries to see if is worth to loop and find the minimum
391 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
392 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
393
394 const Float_t kLimit = 0.02 * 3;
395
396 Float_t dphistarminabs = 1e5;
397 Float_t dphistarmin = 1e5;
398
399 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
400 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
401 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
402 Float_t dphistarabs = TMath::Abs(dphistar);
403
404 if (dphistarabs < dphistarminabs) {
405 dphistarmin = dphistar;
406 dphistarminabs = dphistarabs;
407 }
408 }
409
410 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
411 //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));
412 continue;
413 }
414 }
415 }
416 fHistHBTafter->Fill(deta,dphi);
417 }//HBT cut
418
419 // conversions
420 if(fConversionCut) {
421 if (charge1 * charge2 < 0) {
422 Double_t deta = firstEta - secondEta[j];
423 Double_t dphi = firstPhi - secondPhi[j];
424 fHistConversionbefore->Fill(deta,dphi);
425
426 Float_t m0 = 0.510e-3;
427 Float_t tantheta1 = 1e10;
428
429 // phi in rad
c8ad9635 430 //Float_t phi1rad = firstPhi*TMath::DegToRad();
431 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
432 Float_t phi1rad = firstPhi;
433 Float_t phi2rad = secondPhi[j];
73609a48 434
435 if (firstEta < -1e-10 || firstEta > 1e-10)
436 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
437
438 Float_t tantheta2 = 1e10;
439 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
440 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
441
442 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
443 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
444
445 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
446
447 if (masssqu < 0.04*0.04){
448 //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));
449 continue;
450 }
451 fHistConversionafter->Fill(deta,dphi);
452 }
453 }//conversion cut
9afe3098 454
73609a48 455 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,1.);
456 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,1.);
457 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,1.);
458 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,1.);
8795e993 459 else {
460 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
461 continue;
462 }
6acdbcb2 463 }//end of 2nd particle loop
464 }//end of 1st particle loop
0879e280 465}
466
0879e280 467//____________________________________________________________________//
9afe3098 468TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
469 Int_t iVariablePair,
9afe3098 470 Double_t psiMin,
6acdbcb2 471 Double_t psiMax,
472 Double_t ptTriggerMin,
473 Double_t ptTriggerMax,
474 Double_t ptAssociatedMin,
475 Double_t ptAssociatedMax) {
9afe3098 476 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 477 //(private members) of the AliBalancePsi class.
621329de 478 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
479 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
9afe3098 480 // Psi_2
621329de 481 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
482 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
483 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
484 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
485 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
486 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
9afe3098 487
6acdbcb2 488 // pt trigger
489 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
490 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
491 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
492 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
493 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
494 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
495 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
496 }
497
498 // pt associated
499 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
500 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
501 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
502 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
503 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
504 }
505
a38fd7f3 506 //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));
507
9afe3098 508 // Project into the wanted space (1st: analysis step, 2nd: axis)
509 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
510 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
511 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
512 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
513 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
514 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
515
516 TH1D *gHistBalanceFunctionHistogram = 0x0;
517 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
518 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
519 gHistBalanceFunctionHistogram->Reset();
520
521 switch(iVariablePair) {
621329de 522 case 1:
c8ad9635 523 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
524 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
9afe3098 525 break;
621329de 526 case 2:
c8ad9635 527 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
528 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
9afe3098 529 break;
9afe3098 530 default:
531 break;
532 }
0879e280 533
0879e280 534 hTemp1->Sumw2();
535 hTemp2->Sumw2();
536 hTemp3->Sumw2();
537 hTemp4->Sumw2();
a38fd7f3 538 hTemp1->Add(hTemp3,-1.);
0879e280 539 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 540 hTemp2->Add(hTemp4,-1.);
0879e280 541 hTemp2->Scale(1./hTemp6->GetEntries());
542 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 543 gHistBalanceFunctionHistogram->Scale(0.5);
0879e280 544 }
545
0879e280 546 return gHistBalanceFunctionHistogram;
547}
a38fd7f3 548
f0c5040c 549//____________________________________________________________________//
550TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
551 Int_t iVariablePair,
552 Double_t psiMin,
553 Double_t psiMax,
554 Double_t ptTriggerMin,
555 Double_t ptTriggerMax,
556 Double_t ptAssociatedMin,
557 Double_t ptAssociatedMax,
558 AliBalancePsi *bfMix) {
559 //Returns the BF histogram, extracted from the 6 AliTHn objects
560 //after dividing each correlation function by the Event Mixing one
561 //(private members) of the AliBalancePsi class.
562 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
563 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
564 // Psi_2
565 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
566 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
567 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
568 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
569 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
570 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
571
572 // pt trigger
573 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
574 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
575 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
576 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
577 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
578 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
579 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
580 }
581
582 // pt associated
583 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
584 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
585 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
586 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
587 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
588 }
589
590 // ============================================================================================
591 // the same for event mixing
592 AliTHn *fHistPMix = bfMix->GetHistNp();
593 AliTHn *fHistNMix = bfMix->GetHistNn();
594 AliTHn *fHistPNMix = bfMix->GetHistNpn();
595 AliTHn *fHistNPMix = bfMix->GetHistNnp();
596 AliTHn *fHistPPMix = bfMix->GetHistNpp();
597 AliTHn *fHistNNMix = bfMix->GetHistNnn();
598
599 // Psi_2
600 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
601 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
602 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
603 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
604 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
605 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
606
607 // pt trigger
608 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
609 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
610 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
611 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
612 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
613 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
614 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
615 }
616
617 // pt associated
618 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
619 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
620 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
621 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
622 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
623 }
624 // ============================================================================================
625
626 //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));
627
628 // Project into the wanted space (1st: analysis step, 2nd: axis)
629 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
630 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
631 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
632 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
633 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
634 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
635
636 // ============================================================================================
637 // the same for event mixing
638 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,iVariablePair);
639 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,iVariablePair);
640 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,iVariablePair);
641 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,iVariablePair);
642 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
643 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
644 // ============================================================================================
645
646 TH1D *gHistBalanceFunctionHistogram = 0x0;
647 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
648 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
649 gHistBalanceFunctionHistogram->Reset();
650
651 switch(iVariablePair) {
652 case 1:
653 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
654 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
655 break;
656 case 2:
657 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
658 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
659 break;
660 default:
661 break;
662 }
663
664 hTemp1->Sumw2();
665 hTemp2->Sumw2();
666 hTemp3->Sumw2();
667 hTemp4->Sumw2();
668 hTemp1Mix->Sumw2();
669 hTemp2Mix->Sumw2();
670 hTemp3Mix->Sumw2();
671 hTemp4Mix->Sumw2();
672
673 hTemp1->Scale(1./hTemp5->GetEntries());
674 hTemp3->Scale(1./hTemp5->GetEntries());
675 hTemp2->Scale(1./hTemp6->GetEntries());
676 hTemp4->Scale(1./hTemp6->GetEntries());
677 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
678 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
679 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
680 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
681
682 hTemp1->Divide(hTemp1Mix);
683 hTemp2->Divide(hTemp2Mix);
684 hTemp3->Divide(hTemp3Mix);
685 hTemp4->Divide(hTemp4Mix);
686
687 hTemp1->Add(hTemp3,-1.);
688 hTemp2->Add(hTemp4,-1.);
689
690 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
691 gHistBalanceFunctionHistogram->Scale(0.5);
692 }
693
694 return gHistBalanceFunctionHistogram;
695}
696
a38fd7f3 697//____________________________________________________________________//
621329de 698TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 699 Double_t psiMax,
700 Double_t ptTriggerMin,
701 Double_t ptTriggerMax,
702 Double_t ptAssociatedMin,
703 Double_t ptAssociatedMax) {
621329de 704 //Returns the BF histogram in Delta eta vs Delta phi,
705 //extracted from the 6 AliTHn objects
706 //(private members) of the AliBalancePsi class.
707 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
708 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
709 TString histName = "gHistBalanceFunctionHistogram2D";
710
711 // Psi_2
712 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
713 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
714 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
715 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
716 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
717 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
718
6acdbcb2 719 // pt trigger
720 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
721 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
722 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
723 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
724 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
725 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
726 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
727 }
728
729 // pt associated
730 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
731 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
732 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
733 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
734 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
735 }
736
737 //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 738
739 // Project into the wanted space (1st: analysis step, 2nd: axis)
740 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
741 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
742 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
743 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
744 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
745 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
746
747 TH2D *gHistBalanceFunctionHistogram = 0x0;
748 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
749 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
750 gHistBalanceFunctionHistogram->Reset();
c8ad9635 751 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
752 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
753 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
621329de 754
755 hTemp1->Sumw2();
756 hTemp2->Sumw2();
757 hTemp3->Sumw2();
758 hTemp4->Sumw2();
759 hTemp1->Add(hTemp3,-1.);
760 hTemp1->Scale(1./hTemp5->GetEntries());
761 hTemp2->Add(hTemp4,-1.);
762 hTemp2->Scale(1./hTemp6->GetEntries());
763 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
764 gHistBalanceFunctionHistogram->Scale(0.5);
765 }
766
767 return gHistBalanceFunctionHistogram;
768}
769
f0c5040c 770//____________________________________________________________________//
771TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
772 Double_t psiMax,
773 Double_t ptTriggerMin,
774 Double_t ptTriggerMax,
775 Double_t ptAssociatedMin,
776 Double_t ptAssociatedMax,
777 AliBalancePsi *bfMix) {
778 //Returns the BF histogram in Delta eta vs Delta phi,
779 //after dividing each correlation function by the Event Mixing one
780 //extracted from the 6 AliTHn objects
781 //(private members) of the AliBalancePsi class.
782 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
783 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
784 TString histName = "gHistBalanceFunctionHistogram2D";
785
786 if(!bfMix){
787 AliError("balance function object for event mixing not available");
788 return NULL;
789 }
790
791 // Psi_2
792 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
793 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
794 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
795 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
796 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
797 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
798
799 // pt trigger
800 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
801 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
802 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
803 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
804 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
805 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
806 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
807 }
808
809 // pt associated
810 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
811 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
812 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
813 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
814 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
815 }
816
817
818 // ============================================================================================
819 // the same for event mixing
820 AliTHn *fHistPMix = bfMix->GetHistNp();
821 AliTHn *fHistNMix = bfMix->GetHistNn();
822 AliTHn *fHistPNMix = bfMix->GetHistNpn();
823 AliTHn *fHistNPMix = bfMix->GetHistNnp();
824 AliTHn *fHistPPMix = bfMix->GetHistNpp();
825 AliTHn *fHistNNMix = bfMix->GetHistNnn();
826
827 // Psi_2
828 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
829 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
830 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
831 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
832 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
833 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
834
835 // pt trigger
836 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
837 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
838 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
839 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
840 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
841 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
842 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
843 }
844
845 // pt associated
846 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
847 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
848 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
849 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
850 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
851 }
852 // ============================================================================================
853
854
855 //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)));
856
857 // Project into the wanted space (1st: analysis step, 2nd: axis)
858 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
859 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
860 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
861 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
862 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
863 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
864
865 // ============================================================================================
866 // the same for event mixing
867 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
868 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
869 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
870 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
871 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
872 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
873 // ============================================================================================
874
875 TH2D *gHistBalanceFunctionHistogram = 0x0;
876 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
877 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
878 gHistBalanceFunctionHistogram->Reset();
879 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
880 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
881 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
882
883 hTemp1->Sumw2();
884 hTemp2->Sumw2();
885 hTemp3->Sumw2();
886 hTemp4->Sumw2();
887 hTemp1Mix->Sumw2();
888 hTemp2Mix->Sumw2();
889 hTemp3Mix->Sumw2();
890 hTemp4Mix->Sumw2();
891
892 hTemp1->Scale(1./hTemp5->GetEntries());
893 hTemp3->Scale(1./hTemp5->GetEntries());
894 hTemp2->Scale(1./hTemp6->GetEntries());
895 hTemp4->Scale(1./hTemp6->GetEntries());
896 hTemp1Mix->Scale(1./hTemp5Mix->GetEntries());
897 hTemp3Mix->Scale(1./hTemp5Mix->GetEntries());
898 hTemp2Mix->Scale(1./hTemp6Mix->GetEntries());
899 hTemp4Mix->Scale(1./hTemp6Mix->GetEntries());
900
901 hTemp1->Divide(hTemp1Mix);
902 hTemp2->Divide(hTemp2Mix);
903 hTemp3->Divide(hTemp3Mix);
904 hTemp4->Divide(hTemp4Mix);
905
906 hTemp1->Add(hTemp3,-1.);
907 hTemp2->Add(hTemp4,-1.);
908
909 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
910 gHistBalanceFunctionHistogram->Scale(0.5);
911 }
912
913 return gHistBalanceFunctionHistogram;
914}
915
916
621329de 917//____________________________________________________________________//
918TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 919 Double_t psiMax,
920 Double_t ptTriggerMin,
921 Double_t ptTriggerMax,
922 Double_t ptAssociatedMin,
923 Double_t ptAssociatedMax) {
a38fd7f3 924 //Returns the 2D correlation function for (+-) pairs
621329de 925 // Psi_2: axis 0
926 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
927 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
621329de 928
6acdbcb2 929 // pt trigger
930 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
931 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
932 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
5ef8eaa7 933 }
621329de 934
6acdbcb2 935 // pt associated
936 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
937 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
938
939 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
940 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
941
942 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
943 //TCanvas *c1 = new TCanvas("c1","");
944 //c1->cd();
945 //if(!gHistTest){
946 //AliError("Projection of fHistP = NULL");
947 //return gHistTest;
948 //}
949 //else{
950 //gHistTest->DrawCopy("colz");
951 //}
952
953 //0:step, 1: Delta eta, 2: Delta phi
621329de 954 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 955 if(!gHist){
956 AliError("Projection of fHistPN = NULL");
957 return gHist;
958 }
959
6acdbcb2 960 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
961 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
962 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 963
6acdbcb2 964 //TCanvas *c2 = new TCanvas("c2","");
965 //c2->cd();
966 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 967
968 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
969 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
a38fd7f3 970
971 return gHist;
972}
973
974//____________________________________________________________________//
621329de 975TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 976 Double_t psiMax,
977 Double_t ptTriggerMin,
978 Double_t ptTriggerMax,
979 Double_t ptAssociatedMin,
980 Double_t ptAssociatedMax) {
a38fd7f3 981 //Returns the 2D correlation function for (+-) pairs
621329de 982 // Psi_2: axis 0
983 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
984 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
a38fd7f3 985
6acdbcb2 986 // pt trigger
987 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
988 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
989 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
990 }
991
992 // pt associated
993 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
994 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
995
996 //0:step, 1: Delta eta, 2: Delta phi
621329de 997 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 998 if(!gHist){
999 AliError("Projection of fHistPN = NULL");
1000 return gHist;
1001 }
1002
a38fd7f3 1003 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1004 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
621329de 1005 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1006 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
a38fd7f3 1007
1008 return gHist;
1009}
1010
1011//____________________________________________________________________//
621329de 1012TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 1013 Double_t psiMax,
1014 Double_t ptTriggerMin,
1015 Double_t ptTriggerMax,
1016 Double_t ptAssociatedMin,
1017 Double_t ptAssociatedMax) {
a38fd7f3 1018 //Returns the 2D correlation function for (+-) pairs
621329de 1019 // Psi_2: axis 0
1020 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1021 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
6acdbcb2 1022
1023 // pt trigger
1024 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1025 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1026 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1027 }
1028
1029 // pt associated
1030 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1031 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
a38fd7f3 1032
6acdbcb2 1033 //0:step, 1: Delta eta, 2: Delta phi
621329de 1034 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 1035 if(!gHist){
1036 AliError("Projection of fHistPN = NULL");
1037 return gHist;
1038 }
1039
a38fd7f3 1040 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
1041 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
621329de 1042 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
1043 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
a38fd7f3 1044
1045 return gHist;
1046}
1047
1048//____________________________________________________________________//
621329de 1049TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 1050 Double_t psiMax,
1051 Double_t ptTriggerMin,
1052 Double_t ptTriggerMax,
1053 Double_t ptAssociatedMin,
1054 Double_t ptAssociatedMax) {
a38fd7f3 1055 //Returns the 2D correlation function for (+-) pairs
621329de 1056 // Psi_2: axis 0
1057 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
1058 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
6acdbcb2 1059
1060 // pt trigger
1061 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1062 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1063 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
1064 }
1065
1066 // pt associated
1067 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
1068 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
a38fd7f3 1069
6acdbcb2 1070 //0:step, 1: Delta eta, 2: Delta phi
621329de 1071 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 1072 if(!gHist){
1073 AliError("Projection of fHistPN = NULL");
1074 return gHist;
1075 }
1076
a38fd7f3 1077 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
1078 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
621329de 1079 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
1080 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
a38fd7f3 1081
1082 return gHist;
1083}
621329de 1084
73609a48 1085//____________________________________________________________________//
1086Float_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) {
1087 //
1088 // calculates dphistar
1089 //
1090 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
1091
1092 static const Double_t kPi = TMath::Pi();
1093
1094 // circularity
1095// if (dphistar > 2 * kPi)
1096// dphistar -= 2 * kPi;
1097// if (dphistar < -2 * kPi)
1098// dphistar += 2 * kPi;
1099
1100 if (dphistar > kPi)
1101 dphistar = kPi * 2 - dphistar;
1102 if (dphistar < -kPi)
1103 dphistar = -kPi * 2 - dphistar;
1104 if (dphistar > kPi) // might look funny but is needed
1105 dphistar = kPi * 2 - dphistar;
1106
1107 return dphistar;
1108}
1109