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