]> 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
GetDetTypeName digest both coarse and full det.type
[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),
64 fHistPshiMinusPhi(0),
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),
89 fHistPshiMinusPhi(balance.fHistPshiMinusPhi),
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;
110 delete fHistPshiMinusPhi;
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
6acdbcb2 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 (all)
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;
141 axisTitleSingle[0] = "#phi - #Psi_{2} (a.u.)";
142 iBinPair[0] = kNPsi2Bins;
143 dBinsPair[0] = psi2Bins;
144 axisTitlePair[0] = "#phi - #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;
153 axisTitlePair[1] = "#Delta #eta";
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++){
159 deltaPhiBins[i] = -180.0 + i * 5.;
160 }
621329de 161 iBinPair[2] = kNDeltaPhiBins;
162 dBinsPair[2] = deltaPhiBins;
163 axisTitlePair[2] = "#Delta #phi (#circ)";
9afe3098 164
a38fd7f3 165 // pt(trigger-associated)
8795e993 166 const Int_t kNPtBins = 16;
167 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.};
168 //for(Int_t i = 0; i < kNPtBins+1; i++){
169 //ptBins[i] = 0.2 + i * 0.5;
170 //}
621329de 171 iBinSingle[1] = kNPtBins;
172 dBinsSingle[1] = ptBins;
173 axisTitleSingle[1] = "p_{t}^{trig.} (GeV/c)";
174
175 iBinPair[3] = kNPtBins;
176 dBinsPair[3] = ptBins;
177 axisTitlePair[3] = "p_{t}^{trig.} (GeV/c)";
a38fd7f3 178
179 iBinPair[4] = kNPtBins;
180 dBinsPair[4] = ptBins;
621329de 181 axisTitlePair[4] = "p_{t}^{assoc.} (GeV/c)";
9afe3098 182
183 TString histName;
184 //+ triggered particles
185 histName = "fHistP";
9fd4b54e 186 if(fShuffle) histName.Append("_shuffle");
9afe3098 187 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 188 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
189 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 190 fHistP->SetBinLimits(j, dBinsSingle[j]);
191 fHistP->SetVarTitle(j, axisTitleSingle[j]);
0879e280 192 }
9afe3098 193
194 //- triggered particles
195 histName = "fHistN";
9fd4b54e 196 if(fShuffle) histName.Append("_shuffle");
9afe3098 197 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 198 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
199 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
9afe3098 200 fHistN->SetBinLimits(j, dBinsSingle[j]);
201 fHistN->SetVarTitle(j, axisTitleSingle[j]);
0879e280 202 }
9afe3098 203
204 //+- pairs
205 histName = "fHistPN";
9fd4b54e 206 if(fShuffle) histName.Append("_shuffle");
9afe3098 207 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 208 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
209 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 210 fHistPN->SetBinLimits(j, dBinsPair[j]);
211 fHistPN->SetVarTitle(j, axisTitlePair[j]);
0879e280 212 }
0879e280 213
9afe3098 214 //-+ pairs
215 histName = "fHistNP";
9fd4b54e 216 if(fShuffle) histName.Append("_shuffle");
9afe3098 217 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 218 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
219 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 220 fHistNP->SetBinLimits(j, dBinsPair[j]);
221 fHistNP->SetVarTitle(j, axisTitlePair[j]);
0879e280 222 }
0879e280 223
9afe3098 224 //++ pairs
225 histName = "fHistPP";
9fd4b54e 226 if(fShuffle) histName.Append("_shuffle");
9afe3098 227 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 228 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
229 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 230 fHistPP->SetBinLimits(j, dBinsPair[j]);
231 fHistPP->SetVarTitle(j, axisTitlePair[j]);
232 }
233
234 //-- pairs
235 histName = "fHistNN";
9fd4b54e 236 if(fShuffle) histName.Append("_shuffle");
9afe3098 237 if(fCentralityId) histName += fCentralityId.Data();
9fd4b54e 238 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
239 for (Int_t j=0; j<kTrackVariablesPair; j++) {
9afe3098 240 fHistNN->SetBinLimits(j, dBinsPair[j]);
241 fHistNN->SetVarTitle(j, axisTitlePair[j]);
0879e280 242 }
f06d59b3 243 AliInfo("Finished setting up the AliTHn");
73609a48 244
245 // QA histograms
246 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200);
247 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200);
248 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200);
249 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200);
250 fHistPshiMinusPhi = new TH2D("fHistPshiMinusPhi","",4,-0.5,3.5,100,0,360.);
0879e280 251}
252
253//____________________________________________________________________//
f06d59b3 254void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
255 TObjArray *particles,
73609a48 256 TObjArray *particlesMixed,
257 Float_t bSign) {
0879e280 258 // Calculates the balance function
259 fAnalyzedEvents++;
0879e280 260
261 // Initialize histograms if not done yet
9afe3098 262 if(!fHistPN){
0879e280 263 AliWarning("Histograms not yet initialized! --> Will be done now");
264 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
265 InitHistograms();
266 }
267
9fd4b54e 268 Double_t trackVariablesSingle[kTrackVariablesSingle];
269 Double_t trackVariablesPair[kTrackVariablesPair];
f06d59b3 270
271 if (!particles){
272 AliWarning("particles TObjArray is NULL pointer --> return");
273 return;
274 }
9afe3098 275
f06d59b3 276 // define end of particle loops
277 Int_t iMax = particles->GetEntriesFast();
278 Int_t jMax = iMax;
279 if (particlesMixed)
280 jMax = particlesMixed->GetEntriesFast();
281
282 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
283 TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
284
285 TArrayF secondEta(jMax);
286 TArrayF secondPhi(jMax);
287 TArrayF secondPt(jMax);
288
289 for (Int_t i=0; i<jMax; i++){
290 secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
291 secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
292 secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
293 }
294
295 // 1st particle loop
73609a48 296 for (Int_t i = 0; i < iMax; i++) {
6acdbcb2 297 AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
298
299 // some optimization
300 Float_t firstEta = firstParticle->Eta();
301 Float_t firstPhi = firstParticle->Phi();
302 Float_t firstPt = firstParticle->Pt();
303
304 // Event plane (determine psi bin)
305 Double_t gPsiMinusPhi = 0.;
306 Double_t gPsiMinusPhiBin = -10.;
307 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
308 //in-plane
309 if((gPsiMinusPhi <= 7.5)||
310 ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
311 gPsiMinusPhiBin = 0.0;
312 //intermediate
313 else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
314 ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
315 ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
316 ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
317 gPsiMinusPhiBin = 1.0;
318 //out of plane
319 else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
320 ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
321 gPsiMinusPhiBin = 2.0;
322 //everything else
323 else
324 gPsiMinusPhiBin = 3.0;
325
73609a48 326 fHistPshiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
327
328 Short_t charge1 = (Short_t) firstParticle->Charge();
6acdbcb2 329
330 trackVariablesSingle[0] = gPsiMinusPhiBin;
331 trackVariablesSingle[1] = firstPt;
332
333 //fill single particle histograms
73609a48 334 if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,1.);
335 else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,1.);
6acdbcb2 336
337 // 2nd particle loop (only for j < i for non double counting in the same pT region)
338 // --> SAME pT region for trigger and assoc: NO double counting with this
339 // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i
340 // --> can be handled afterwards by using assoc. as trigger as well ?!
a86235ae 341 for(Int_t j = 0; j < iMax; j++) {
342
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);
9afe3098 349
6acdbcb2 350 Short_t charge2 = (Short_t) secondParticle->Charge();
9afe3098 351
6acdbcb2 352 trackVariablesPair[0] = gPsiMinusPhiBin;
353 trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
354 trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
355 if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
356 trackVariablesPair[2] -= 360.;
357 if (trackVariablesPair[2] < - 180.)
358 trackVariablesPair[2] += 360.;
9afe3098 359
6acdbcb2 360 trackVariablesPair[3] = firstPt; // pt trigger
361 trackVariablesPair[4] = secondPt[j]; // pt
362 // trackVariablesPair[5] = fCentrality; // centrality
73609a48 363
364 // HBT like cut
365 if(fHBTCut && charge1 * charge2 > 0){
366 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
367 // continue;
368
369 Double_t deta = firstEta - secondEta[j];
370 Double_t dphi = firstPhi - secondPhi[j];
371 // VERSION 2 (Taken from DPhiCorrelations)
372 // the variables & cuthave been developed by the HBT group
373 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
374 fHistHBTbefore->Fill(deta,dphi);
375
376 // optimization
377 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
378 {
379 // phi in rad
380 Float_t phi1rad = firstPhi*TMath::DegToRad();
381 Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
382
383 // check first boundaries to see if is worth to loop and find the minimum
384 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
385 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
386
387 const Float_t kLimit = 0.02 * 3;
388
389 Float_t dphistarminabs = 1e5;
390 Float_t dphistarmin = 1e5;
391
392 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
393 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
394 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
395 Float_t dphistarabs = TMath::Abs(dphistar);
396
397 if (dphistarabs < dphistarminabs) {
398 dphistarmin = dphistar;
399 dphistarminabs = dphistarabs;
400 }
401 }
402
403 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02) {
404 //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));
405 continue;
406 }
407 }
408 }
409 fHistHBTafter->Fill(deta,dphi);
410 }//HBT cut
411
412 // conversions
413 if(fConversionCut) {
414 if (charge1 * charge2 < 0) {
415 Double_t deta = firstEta - secondEta[j];
416 Double_t dphi = firstPhi - secondPhi[j];
417 fHistConversionbefore->Fill(deta,dphi);
418
419 Float_t m0 = 0.510e-3;
420 Float_t tantheta1 = 1e10;
421
422 // phi in rad
423 Float_t phi1rad = firstPhi*TMath::DegToRad();
424 Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
425
426 if (firstEta < -1e-10 || firstEta > 1e-10)
427 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
428
429 Float_t tantheta2 = 1e10;
430 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
431 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
432
433 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
434 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
435
436 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
437
438 if (masssqu < 0.04*0.04){
439 //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));
440 continue;
441 }
442 fHistConversionafter->Fill(deta,dphi);
443 }
444 }//conversion cut
9afe3098 445
73609a48 446 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,1.);
447 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,1.);
448 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,1.);
449 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,1.);
8795e993 450 else {
451 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
452 continue;
453 }
6acdbcb2 454 }//end of 2nd particle loop
455 }//end of 1st particle loop
0879e280 456}
457
0879e280 458//____________________________________________________________________//
9afe3098 459TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
460 Int_t iVariablePair,
9afe3098 461 Double_t psiMin,
6acdbcb2 462 Double_t psiMax,
463 Double_t ptTriggerMin,
464 Double_t ptTriggerMax,
465 Double_t ptAssociatedMin,
466 Double_t ptAssociatedMax) {
9afe3098 467 //Returns the BF histogram, extracted from the 6 AliTHn objects
0879e280 468 //(private members) of the AliBalancePsi class.
621329de 469 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
470 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
9afe3098 471 // Psi_2
621329de 472 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
473 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
474 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
475 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
476 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
477 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
9afe3098 478
6acdbcb2 479 // pt trigger
480 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
481 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
482 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
483 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
484 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
485 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
486 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
487 }
488
489 // pt associated
490 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
491 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
492 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
493 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
494 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
495 }
496
a38fd7f3 497 //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));
498
9afe3098 499 // Project into the wanted space (1st: analysis step, 2nd: axis)
500 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair);
501 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair);
502 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair);
503 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair);
504 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
505 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
506
507 TH1D *gHistBalanceFunctionHistogram = 0x0;
508 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
509 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
510 gHistBalanceFunctionHistogram->Reset();
511
512 switch(iVariablePair) {
621329de 513 case 1:
9afe3098 514 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
515 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
516 break;
621329de 517 case 2:
9afe3098 518 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
519 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
520 break;
9afe3098 521 default:
522 break;
523 }
0879e280 524
0879e280 525 hTemp1->Sumw2();
526 hTemp2->Sumw2();
527 hTemp3->Sumw2();
528 hTemp4->Sumw2();
a38fd7f3 529 hTemp1->Add(hTemp3,-1.);
0879e280 530 hTemp1->Scale(1./hTemp5->GetEntries());
a38fd7f3 531 hTemp2->Add(hTemp4,-1.);
0879e280 532 hTemp2->Scale(1./hTemp6->GetEntries());
533 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
9afe3098 534 gHistBalanceFunctionHistogram->Scale(0.5);
0879e280 535 }
536
0879e280 537 return gHistBalanceFunctionHistogram;
538}
a38fd7f3 539
540//____________________________________________________________________//
621329de 541TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
6acdbcb2 542 Double_t psiMax,
543 Double_t ptTriggerMin,
544 Double_t ptTriggerMax,
545 Double_t ptAssociatedMin,
546 Double_t ptAssociatedMax) {
621329de 547 //Returns the BF histogram in Delta eta vs Delta phi,
548 //extracted from the 6 AliTHn objects
549 //(private members) of the AliBalancePsi class.
550 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
551 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
552 TString histName = "gHistBalanceFunctionHistogram2D";
553
554 // Psi_2
555 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
556 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
557 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
558 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
559 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
560 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
561
6acdbcb2 562 // pt trigger
563 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
564 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
565 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
566 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
567 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
568 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
569 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
570 }
571
572 // pt associated
573 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
574 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
575 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
576 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
577 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
578 }
579
580 //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 581
582 // Project into the wanted space (1st: analysis step, 2nd: axis)
583 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
584 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
585 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
586 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
587 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
588 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
589
590 TH2D *gHistBalanceFunctionHistogram = 0x0;
591 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
592 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
593 gHistBalanceFunctionHistogram->Reset();
594 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
595 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta #phi (deg.)");
596 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta #eta,#Delta #phi)");
597
598 hTemp1->Sumw2();
599 hTemp2->Sumw2();
600 hTemp3->Sumw2();
601 hTemp4->Sumw2();
602 hTemp1->Add(hTemp3,-1.);
603 hTemp1->Scale(1./hTemp5->GetEntries());
604 hTemp2->Add(hTemp4,-1.);
605 hTemp2->Scale(1./hTemp6->GetEntries());
606 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
607 gHistBalanceFunctionHistogram->Scale(0.5);
608 }
609
610 return gHistBalanceFunctionHistogram;
611}
612
613//____________________________________________________________________//
614TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
6acdbcb2 615 Double_t psiMax,
616 Double_t ptTriggerMin,
617 Double_t ptTriggerMax,
618 Double_t ptAssociatedMin,
619 Double_t ptAssociatedMax) {
a38fd7f3 620 //Returns the 2D correlation function for (+-) pairs
621329de 621 // Psi_2: axis 0
622 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
623 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
621329de 624
6acdbcb2 625 // pt trigger
626 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
627 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
628 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
5ef8eaa7 629 }
621329de 630
6acdbcb2 631 // pt associated
632 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
633 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
634
635 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
636 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
637
638 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
639 //TCanvas *c1 = new TCanvas("c1","");
640 //c1->cd();
641 //if(!gHistTest){
642 //AliError("Projection of fHistP = NULL");
643 //return gHistTest;
644 //}
645 //else{
646 //gHistTest->DrawCopy("colz");
647 //}
648
649 //0:step, 1: Delta eta, 2: Delta phi
621329de 650 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
89c00c43 651 if(!gHist){
652 AliError("Projection of fHistPN = NULL");
653 return gHist;
654 }
655
6acdbcb2 656 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
657 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
658 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
621329de 659
6acdbcb2 660 //TCanvas *c2 = new TCanvas("c2","");
661 //c2->cd();
662 //fHistPN->Project(0,1,2)->DrawCopy("colz");
621329de 663
664 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
665 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
a38fd7f3 666
667 return gHist;
668}
669
670//____________________________________________________________________//
621329de 671TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
6acdbcb2 672 Double_t psiMax,
673 Double_t ptTriggerMin,
674 Double_t ptTriggerMax,
675 Double_t ptAssociatedMin,
676 Double_t ptAssociatedMax) {
a38fd7f3 677 //Returns the 2D correlation function for (+-) pairs
621329de 678 // Psi_2: axis 0
679 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
680 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
a38fd7f3 681
6acdbcb2 682 // pt trigger
683 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
684 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
685 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
686 }
687
688 // pt associated
689 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
690 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
691
692 //0:step, 1: Delta eta, 2: Delta phi
621329de 693 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
89c00c43 694 if(!gHist){
695 AliError("Projection of fHistPN = NULL");
696 return gHist;
697 }
698
a38fd7f3 699 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
700 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
621329de 701 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
702 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
a38fd7f3 703
704 return gHist;
705}
706
707//____________________________________________________________________//
621329de 708TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
6acdbcb2 709 Double_t psiMax,
710 Double_t ptTriggerMin,
711 Double_t ptTriggerMax,
712 Double_t ptAssociatedMin,
713 Double_t ptAssociatedMax) {
a38fd7f3 714 //Returns the 2D correlation function for (+-) pairs
621329de 715 // Psi_2: axis 0
716 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
717 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
6acdbcb2 718
719 // pt trigger
720 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
721 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
722 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
723 }
724
725 // pt associated
726 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
727 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
a38fd7f3 728
6acdbcb2 729 //0:step, 1: Delta eta, 2: Delta phi
621329de 730 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
89c00c43 731 if(!gHist){
732 AliError("Projection of fHistPN = NULL");
733 return gHist;
734 }
735
a38fd7f3 736 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
737 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
621329de 738 if((Double_t)(fHistP->Project(0,1)->GetEntries())!=0)
739 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->GetEntries()));
a38fd7f3 740
741 return gHist;
742}
743
744//____________________________________________________________________//
621329de 745TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
6acdbcb2 746 Double_t psiMax,
747 Double_t ptTriggerMin,
748 Double_t ptTriggerMax,
749 Double_t ptAssociatedMin,
750 Double_t ptAssociatedMax) {
a38fd7f3 751 //Returns the 2D correlation function for (+-) pairs
621329de 752 // Psi_2: axis 0
753 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
754 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax);
6acdbcb2 755
756 // pt trigger
757 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
758 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax);
759 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax);
760 }
761
762 // pt associated
763 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
764 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax);
a38fd7f3 765
6acdbcb2 766 //0:step, 1: Delta eta, 2: Delta phi
621329de 767 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
89c00c43 768 if(!gHist){
769 AliError("Projection of fHistPN = NULL");
770 return gHist;
771 }
772
a38fd7f3 773 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
774 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
621329de 775 if((Double_t)(fHistN->Project(0,1)->GetEntries())!=0)
776 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->GetEntries()));
a38fd7f3 777
778 return gHist;
779}
621329de 780
73609a48 781//____________________________________________________________________//
782Float_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) {
783 //
784 // calculates dphistar
785 //
786 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
787
788 static const Double_t kPi = TMath::Pi();
789
790 // circularity
791// if (dphistar > 2 * kPi)
792// dphistar -= 2 * kPi;
793// if (dphistar < -2 * kPi)
794// dphistar += 2 * kPi;
795
796 if (dphistar > kPi)
797 dphistar = kPi * 2 - dphistar;
798 if (dphistar < -kPi)
799 dphistar = -kPi * 2 - dphistar;
800 if (dphistar > kPi) // might look funny but is needed
801 dphistar = kPi * 2 - dphistar;
802
803 return dphistar;
804}
805