]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/BalanceFunctions/AliBalance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
activate wider tpc cut
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalance.cxx
CommitLineData
6c178944 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$ */
15
16//-----------------------------------------------------------------
17// Balance Function class
18// This is the class to deal with the Balance Function analysis
19// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
20//-----------------------------------------------------------------
21
22
23//ROOT
24#include <Riostream.h>
25#include <TMath.h>
d560b581 26#include <TAxis.h>
74b04dd4 27#include <TH2D.h>
6c178944 28#include <TLorentzVector.h>
5c33329d 29#include <TObjArray.h>
6c178944 30#include <TGraphErrors.h>
9d1f0df5 31#include <TString.h>
6c178944 32
cd54a838 33#include "AliVParticle.h"
34#include "AliMCParticle.h"
5c33329d 35#include "AliESDtrack.h"
cd54a838 36#include "AliAODTrack.h"
37
6c178944 38#include "AliBalance.h"
39
40ClassImp(AliBalance)
41
7f0257ea 42//____________________________________________________________________//
5d534bd4 43AliBalance::AliBalance() :
44 TObject(),
6432ac6a 45 bShuffle(kFALSE),
49a3750f 46 bHBTcut(kFALSE),
845fdeca 47 bConversionCut(kFALSE),
3feee083 48 fAnalysisLevel("ESD"),
6a3f819e 49 fAnalyzedEvents(0) ,
74b04dd4 50 fCentralityId(0) ,
51 fCentStart(0.),
52 fCentStop(0.)
53{
6c178944 54 // Default constructor
3feee083 55
b966652f 56 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
57 if(i == 6) {
58 fNumberOfBins[i] = 180;
6432ac6a 59 fP1Start[i] = -360.0;
b966652f 60 fP1Stop[i] = 360.0;
61 fP2Start[i] = -360.0;
62 fP2Stop[i] = 360.0;
63 fP2Step[i] = 0.1;
64 }
65 else {
66 fNumberOfBins[i] = 20;
67 fP1Start[i] = -1.0;
68 fP1Stop[i] = 1.0;
69 fP2Start[i] = 0.0;
70 fP2Stop[i] = 2.0;
71 }
72 fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
74b04dd4 73 fCentStart = 0.;
74 fCentStop = 0.;
75
b966652f 76 fNn[i] = 0.0;
77 fNp[i] = 0.0;
3feee083 78
79 for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
80 fNpp[i][j] = .0;
81 fNnn[i][j] = .0;
82 fNpn[i][j] = .0;
6432ac6a 83 fNnp[i][j] = .0;
3feee083 84 fB[i][j] = 0.0;
85 ferror[i][j] = 0.0;
6432ac6a 86 }
b966652f 87
6432ac6a 88 fHistP[i] = NULL;
89 fHistN[i] = NULL;
90 fHistPP[i] = NULL;
91 fHistPN[i] = NULL;
92 fHistNP[i] = NULL;
93 fHistNN[i] = NULL;
b966652f 94
95 }
6c178944 96}
97
6c178944 98
7f0257ea 99//____________________________________________________________________//
100AliBalance::AliBalance(const AliBalance& balance):
49a3750f 101 TObject(balance),
102 bShuffle(balance.bShuffle),
103 bHBTcut(balance.bHBTcut),
845fdeca 104 bConversionCut(balance.bConversionCut),
6a3f819e 105 fAnalysisLevel(balance.fAnalysisLevel),
106 fAnalyzedEvents(balance.fAnalyzedEvents),
74b04dd4 107 fCentralityId(balance.fCentralityId),
108 fCentStart(balance.fCentStart),
109 fCentStop(balance.fCentStop) {
7f0257ea 110 //copy constructor
b966652f 111 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
112 fNn[i] = balance.fNn[i];
113 fNp[i] = balance.fNp[i];
3feee083 114
b966652f 115 fP1Start[i] = balance.fP1Start[i];
116 fP1Stop[i] = balance.fP1Stop[i];
3feee083 117 fNumberOfBins[i] = balance.fNumberOfBins[i];
118 fP2Start[i] = balance.fP2Start[i];
119 fP2Stop[i] = balance.fP2Stop[i];
120 fP2Step[i] = balance.fP2Step[i];
74b04dd4 121 fCentStart = balance.fCentStart;
122 fCentStop = balance.fCentStop;
3feee083 123
b966652f 124 fHistP[i] = balance.fHistP[i];
125 fHistN[i] = balance.fHistN[i];
126 fHistPN[i] = balance.fHistPN[i];
127 fHistNP[i] = balance.fHistNP[i];
128 fHistPP[i] = balance.fHistPP[i];
129 fHistNN[i] = balance.fHistNN[i];
130
3feee083 131 for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
132 fNpp[i][j] = .0;
133 fNnn[i][j] = .0;
134 fNpn[i][j] = .0;
6432ac6a 135 fNnp[i][j] = .0;
3feee083 136 fB[i][j] = 0.0;
137 ferror[i][j] = 0.0;
138 }
139 }
140 }
141
7f0257ea 142
143//____________________________________________________________________//
144AliBalance::~AliBalance() {
6c178944 145 // Destructor
6c178944 146
6432ac6a 147 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
148
149 delete fHistP[i];
150 delete fHistN[i];
151 delete fHistPN[i];
152 delete fHistNP[i];
153 delete fHistPP[i];
154 delete fHistNN[i];
155
156 }
7f0257ea 157}
158
7f0257ea 159//____________________________________________________________________//
6432ac6a 160void AliBalance::SetInterval(Int_t iAnalysisType,
161 Double_t p1Start, Double_t p1Stop,
162 Int_t ibins, Double_t p2Start, Double_t p2Stop) {
6c178944 163 // Sets the analyzed interval.
3feee083 164 // Set the same Information for all analyses
6432ac6a 165
166 if(iAnalysisType == -1){
b966652f 167 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
168 fP1Start[i] = p1Start;
169 fP1Stop[i] = p1Stop;
6432ac6a 170 fNumberOfBins[i] = ibins;
3feee083 171 fP2Start[i] = p2Start;
172 fP2Stop[i] = p2Stop;
173 fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
174 }
7f0257ea 175 }
3feee083 176 // Set the Information for one analysis
6432ac6a 177 else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
178 fP1Start[iAnalysisType] = p1Start;
179 fP1Stop[iAnalysisType] = p1Stop;
180 fNumberOfBins[iAnalysisType] = ibins;
181 fP2Start[iAnalysisType] = p2Start;
182 fP2Stop[iAnalysisType] = p2Stop;
183 fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
7f0257ea 184 }
6432ac6a 185 else {
3feee083 186 AliError("Wrong ANALYSIS number!");
7f0257ea 187 }
6c178944 188}
189
6432ac6a 190//____________________________________________________________________//
191void AliBalance::InitHistograms() {
6a3f819e 192 //Initialize the histograms
6432ac6a 193 TString histName;
194 for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
6a3f819e 195 histName = "fHistP"; histName += gBFAnalysisType[iAnalysisType];
196 if(bShuffle) histName.Append("_shuffle");
197 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 198 fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
6a3f819e 199
200 histName = "fHistN"; histName += gBFAnalysisType[iAnalysisType];
201 if(bShuffle) histName.Append("_shuffle");
202 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 203 fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
6432ac6a 204
6a3f819e 205 histName = "fHistPN"; histName += gBFAnalysisType[iAnalysisType];
206 if(bShuffle) histName.Append("_shuffle");
207 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 208 fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
6a3f819e 209
210 histName = "fHistNP"; histName += gBFAnalysisType[iAnalysisType];
211 if(bShuffle) histName.Append("_shuffle");
212 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 213 fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
6a3f819e 214
215 histName = "fHistPP"; histName += gBFAnalysisType[iAnalysisType];
216 if(bShuffle) histName.Append("_shuffle");
217 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 218 fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
6a3f819e 219
220 histName = "fHistNN"; histName += gBFAnalysisType[iAnalysisType];
221 if(bShuffle) histName.Append("_shuffle");
222 if(fCentralityId) histName += fCentralityId.Data();
74b04dd4 223 fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
6432ac6a 224 }
225}
3feee083 226
898d5b28 227//____________________________________________________________________//
09bb7bf4 228void AliBalance::PrintAnalysisSettings() {
09bb7bf4 229
230 Printf("======================================");
231 Printf("Analysis level: %s",fAnalysisLevel.Data());
3feee083 232 Printf("======================================");
b966652f 233 for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
3feee083 234 Printf("Interval info for variable %d",ibin);
235 Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
236 Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
237 Printf("Number of bins: %d",fNumberOfBins[ibin]);
238 Printf("Step: %lf",fP2Step[ibin]);
239 Printf(" ");
240 }
09bb7bf4 241 Printf("======================================");
898d5b28 242}
243
7f0257ea 244//____________________________________________________________________//
845fdeca 245void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Float_t bSign) {
6c178944 246 // Calculates the balance function
247 fAnalyzedEvents++;
248 Int_t i = 0 , j = 0;
3feee083 249 Int_t iBin = 0;
845fdeca 250
6432ac6a 251 // Initialize histograms if not done yet
252 if(!fHistPN[0]){
253 AliWarning("Histograms not yet initialized! --> Will be done now");
254 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
255 InitHistograms();
256 }
257
74b04dd4 258 Int_t gNtrack = chargeVector[0]->size();
259 //Printf("(AliBalance) Number of tracks: %d",gNtrack);
260
261 for(i = 0; i < gNtrack;i++){
f12385bb 262 Short_t charge = chargeVector[0]->at(i);
74b04dd4 263 Double_t rapidity = chargeVector[1]->at(i);
264 Double_t pseudorapidity = chargeVector[2]->at(i);
265 Double_t phi = chargeVector[3]->at(i);
b966652f 266
267 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
268 for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
6432ac6a 269 if(iAnalysisType == kEta) {
b966652f 270 if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
271 if(charge > 0) {
272 fNp[iAnalysisType] += 1.;
74b04dd4 273 fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
b966652f 274 }//charge > 0
275 if(charge < 0) {
276 fNn[iAnalysisType] += 1.;
74b04dd4 277 fHistN[iAnalysisType]->Fill(fCentrality,pseudorapidity);
b966652f 278 }//charge < 0
279 }//p1 interval check
280 }//analysis type: eta
5d42ac40 281 else if(iAnalysisType == kPhi) {
b966652f 282 if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
283 if(charge > 0) {
284 fNp[iAnalysisType] += 1.;
74b04dd4 285 fHistP[iAnalysisType]->Fill(fCentrality,phi);
b966652f 286 }//charge > 0
287 if(charge < 0) {
288 fNn[iAnalysisType] += 1.;
74b04dd4 289 fHistN[iAnalysisType]->Fill(fCentrality,phi);
b966652f 290 }//charge < 0
291 }//p1 interval check
292 }//analysis type: phi
293 else {
294 if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
295 if(charge > 0) {
296 fNp[iAnalysisType] += 1.;
74b04dd4 297 fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
b966652f 298 }//charge > 0
299 if(charge < 0) {
300 fNn[iAnalysisType] += 1.;
74b04dd4 301 fHistN[iAnalysisType]->Fill(fCentrality,rapidity);
b966652f 302 }//charge < 0
303 }//p1 interval check
304 }//analysis type: y, qside, qout, qlong, qinv
305 }//analysis type loop
7f0257ea 306 }
74b04dd4 307
32a4f715 308 //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);
6c178944 309
b966652f 310 Double_t dy = 0., deta = 0.;
311 Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
312 Double_t dphi = 0.;
313
f12385bb 314 Short_t charge1 = 0;
b966652f 315 Double_t eta1 = 0., rap1 = 0.;
74b04dd4 316 Double_t px1 = 0., py1 = 0., pz1 = 0.;
317 Double_t pt1 = 0.;
b966652f 318 Double_t energy1 = 0.;
319 Double_t phi1 = 0.;
320
f12385bb 321 Short_t charge2 = 0;
b966652f 322 Double_t eta2 = 0., rap2 = 0.;
74b04dd4 323 Double_t px2 = 0., py2 = 0., pz2 = 0.;
324 Double_t pt2 = 0.;
b966652f 325 Double_t energy2 = 0.;
326 Double_t phi2 = 0.;
6c178944 327 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
3feee083 328 for(i = 1; i < gNtrack; i++) {
c511520f 329
74b04dd4 330 charge1 = chargeVector[0]->at(i);
331 rap1 = chargeVector[1]->at(i);
332 eta1 = chargeVector[2]->at(i);
333 phi1 = chargeVector[3]->at(i);
334 px1 = chargeVector[4]->at(i);
335 py1 = chargeVector[5]->at(i);
336 pz1 = chargeVector[6]->at(i);
337 pt1 = chargeVector[7]->at(i);
338 energy1 = chargeVector[8]->at(i);
3feee083 339
340 for(j = 0; j < i; j++) {
74b04dd4 341
342 charge2 = chargeVector[0]->at(j);
343 rap2 = chargeVector[1]->at(j);
344 eta2 = chargeVector[2]->at(j);
345 phi2 = chargeVector[3]->at(j);
346 px2 = chargeVector[4]->at(j);
347 py2 = chargeVector[5]->at(j);
348 pz2 = chargeVector[6]->at(j);
349 pt2 = chargeVector[7]->at(i);
350 energy2 = chargeVector[8]->at(j);
351
3feee083 352 // filling the arrays
353
354 // RAPIDITY
74b04dd4 355 dy = TMath::Abs(rap1 - rap2);
b966652f 356
3feee083 357 // Eta
b966652f 358 deta = TMath::Abs(eta1 - eta2);
3feee083 359
b966652f 360 //qlong
3feee083 361 Double_t eTot = energy1 + energy2;
74b04dd4 362 Double_t pxTot = px1 + px2;
363 Double_t pyTot = py1 + py2;
364 Double_t pzTot = pz1 + pz2;
3feee083 365 Double_t q0Tot = energy1 - energy2;
74b04dd4 366 Double_t qxTot = px1 - px2;
367 Double_t qyTot = py1 - py2;
368 Double_t qzTot = pz1 - pz2;
369
370 Double_t eTot2 = eTot*eTot;
371 Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
372 Double_t pzTot2 = pzTot*pzTot;
373
374 Double_t q0Tot2 = q0Tot*q0Tot;
375 Double_t qTot2 = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
376
377 Double_t snn = eTot2 - pTot2;
378 Double_t ptTot2 = pTot2 - pzTot2 ;
379 Double_t ptTot = TMath::Sqrt( ptTot2 );
3feee083 380
74b04dd4 381 qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
b966652f 382
383 //qout
74b04dd4 384 qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
b966652f 385
386 //qside
387 qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
3feee083 388
b966652f 389 //qinv
74b04dd4 390 qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
b966652f 391
392 //phi
6432ac6a 393 dphi = TMath::Abs(phi1 - phi2);
4bd2c3ed 394 if(dphi>180) dphi = 360 - dphi; //dphi should be between 0 and 180!
b966652f 395
49a3750f 396 // HBT like cut
397 if(bHBTcut){
845fdeca 398 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
399 // continue;
400
401 // VERSION 2 (Taken from DPhiCorrelations)
402 // the variables & cuthave been developed by the HBT group
403 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
404
405 // optimization
406 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
407 {
408 // check first boundaries to see if is worth to loop and find the minimum
409 Float_t dphistar1 = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, 0.8, bSign);
410 Float_t dphistar2 = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, 2.5, bSign);
411
412 const Float_t kLimit = 0.02 * 3;
413
414 Float_t dphistarminabs = 1e5;
415 Float_t dphistarmin = 1e5;
416 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
417 {
418 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
419 {
420 Float_t dphistar = GetDPhiStar(phi1*TMath::DegToRad(), pt1, charge1, phi2*TMath::DegToRad(), pt2, charge2, rad, bSign);
421
422 Float_t dphistarabs = TMath::Abs(dphistar);
423
424 if (dphistarabs < dphistarminabs)
425 {
426 dphistarmin = dphistar;
427 dphistarminabs = dphistarabs;
428 }
429 }
430
431 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
432 {
433 //Printf("Removed track pair %d %d with %f %f %f %f %d %f %f %d %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
434 continue;
435 }
436 }
437 }
49a3750f 438 }
845fdeca 439
440 // conversions
441 if(bConversionCut){
442 if (charge1 * charge2 < 0)
443 {
444 Float_t m0 = 0.510e-3;
445 Float_t tantheta1 = 1e10;
446
447 if (eta1 < -1e-10 || eta1 > 1e-10)
448 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
449
450 Float_t tantheta2 = 1e10;
451 if (eta2 < -1e-10 || eta2 > 1e-10)
452 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
453
454 Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
455 Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
456
457 Float_t mass = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
458
459 if (mass < 0.04*0.04){
460 //Printf("Removed track pair %d %d with %f %f %d %d ", i, j, deta, mass, charge1, charge2);
461 continue;
462 }
463 }
464 }
465
466
b966652f 467 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
74b04dd4 468 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
469
470 // rapidity
471 if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
472 iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
473 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
474
475 if((charge1 > 0)&&(charge2 > 0)) {
476 fNpp[kRapidity][iBin] += 1.;
477 fHistPP[kRapidity]->Fill(fCentrality,dy);
478 }
479 else if((charge1 < 0)&&(charge2 < 0)) {
480 fNnn[kRapidity][iBin] += 1.;
481 fHistNN[kRapidity]->Fill(fCentrality,dy);
482 }
483 else if((charge1 > 0)&&(charge2 < 0)) {
484 fNpn[kRapidity][iBin] += 1.;
485 fHistPN[kRapidity]->Fill(fCentrality,dy);
486 }
487 else if((charge1 < 0)&&(charge2 > 0)) {
488 fNpn[kRapidity][iBin] += 1.;
489 fHistPN[kRapidity]->Fill(fCentrality,dy);
490 }
491 }//BF binning check
492 }//p2 interval check
493 }//p1 interval check
6432ac6a 494
74b04dd4 495 // pseudorapidity
496 if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
497 if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
498 iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
499 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
500 if((charge1 > 0)&&(charge2 > 0)) {
501 fNpp[kEta][iBin] += 1.;
502 fHistPP[kEta]->Fill(fCentrality,deta);
503 }
504 if((charge1 < 0)&&(charge2 < 0)) {
505 fNnn[kEta][iBin] += 1.;
506 fHistNN[kEta]->Fill(fCentrality,deta);
507 }
508 if((charge1 > 0)&&(charge2 < 0)) {
509 fNpn[kEta][iBin] += 1.;
510 fHistPN[kEta]->Fill(fCentrality,deta);
511 }
512 if((charge1 < 0)&&(charge2 > 0)) {
513 fNpn[kEta][iBin] += 1.;
514 fHistPN[kEta]->Fill(fCentrality,deta);
515 }
516 }//BF binning check
517 }//p2 interval check
518 }//p1 interval check
b966652f 519
74b04dd4 520 // Qlong, out, side, inv
521 // Check the p1 intervall for rapidity here (like for single tracks above)
522 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
523 if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
524 iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
525 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
526 if((charge1 > 0)&&(charge2 > 0)) {
527 fNpp[kQlong][iBin] += 1.;
528 fHistPP[kQlong]->Fill(fCentrality,qLong);
529 }
530 if((charge1 < 0)&&(charge2 < 0)) {
531 fNnn[kQlong][iBin] += 1.;
532 fHistNN[kQlong]->Fill(fCentrality,qLong);
533 }
534 if((charge1 > 0)&&(charge2 < 0)) {
535 fNpn[kQlong][iBin] += 1.;
536 fHistPN[kQlong]->Fill(fCentrality,qLong);
537 }
538 if((charge1 < 0)&&(charge2 > 0)) {
539 fNpn[kQlong][iBin] += 1.;
540 fHistPN[kQlong]->Fill(fCentrality,qLong);
541 }
542 }//BF binning check
543 }//p2 interval check
544 }//p1 interval check
6432ac6a 545
74b04dd4 546 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
547 if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
548 iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
549 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
550 if((charge1 > 0)&&(charge2 > 0)) {
551 fNpp[kQout][iBin] += 1.;
552 fHistPP[kQout]->Fill(fCentrality,qOut);
553 }
554 if((charge1 < 0)&&(charge2 < 0)) {
555 fNnn[kQout][iBin] += 1.;
556 fHistNN[kQout]->Fill(fCentrality,qOut);
557 }
558 if((charge1 > 0)&&(charge2 < 0)) {
559 fNpn[kQout][iBin] += 1.;
560 fHistPN[kQout]->Fill(fCentrality,qOut);
561 }
562 if((charge1 < 0)&&(charge2 > 0)) {
563 fNpn[kQout][iBin] += 1.;
564 fHistPN[kQout]->Fill(fCentrality,qOut);
565 }
566 }//BF binning check
567 }//p2 interval check
568 }//p1 interval check
569
570 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
571 if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
572 iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
573 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
574 if((charge1 > 0)&&(charge2 > 0)) {
575 fNpp[kQside][iBin] += 1.;
576 fHistPP[kQside]->Fill(fCentrality,qSide);
577 }
578 if((charge1 < 0)&&(charge2 < 0)) {
579 fNnn[kQside][iBin] += 1.;
580 fHistNN[kQside]->Fill(fCentrality,qSide);
581 }
582 if((charge1 > 0)&&(charge2 < 0)) {
583 fNpn[kQside][iBin] += 1.;
584 fHistPN[kQside]->Fill(fCentrality,qSide);
585 }
586 if((charge1 < 0)&&(charge2 > 0)) {
587 fNpn[kQside][iBin] += 1.;
588 fHistPN[kQside]->Fill(fCentrality,qSide);
589 }
590 }//BF binning check
591 }//p2 interval check
592 }//p1 interval check
593
594 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
595 if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
596 iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
597 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
598 if((charge1 > 0)&&(charge2 > 0)) {
599 fNpp[kQinv][iBin] += 1.;
600 fHistPP[kQinv]->Fill(fCentrality,qInv);
601 }
602 if((charge1 < 0)&&(charge2 < 0)) {
603 fNnn[kQinv][iBin] += 1.;
604 fHistNN[kQinv]->Fill(fCentrality,qInv);
605 }
606 if((charge1 > 0)&&(charge2 < 0)) {
607 fNpn[kQinv][iBin] += 1.;
608 fHistPN[kQinv]->Fill(fCentrality,qInv);
609 }
610 if((charge1 < 0)&&(charge2 > 0)) {
611 fNpn[kQinv][iBin] += 1.;
612 fHistPN[kQinv]->Fill(fCentrality,qInv);
613 }
614 }//BF binning check
615 }//p2 interval check
616 }//p1 interval check
617
618 // Phi
619 if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
620 if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
621 iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
622 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
623 if((charge1 > 0)&&(charge2 > 0)) {
624 fNpp[kPhi][iBin] += 1.;
625 fHistPP[kPhi]->Fill(fCentrality,dphi);
626 }
627 if((charge1 < 0)&&(charge2 < 0)) {
628 fNnn[kPhi][iBin] += 1.;
629 fHistNN[kPhi]->Fill(fCentrality,dphi);
630 }
631 if((charge1 > 0)&&(charge2 < 0)) {
632 fNpn[kPhi][iBin] += 1.;
633 fHistPN[kPhi]->Fill(fCentrality,dphi);
634 }
635 if((charge1 < 0)&&(charge2 > 0)) {
636 fNpn[kPhi][iBin] += 1.;
637 fHistPN[kPhi]->Fill(fCentrality,dphi);
638 }
639 }//BF binning check
640 }//p2 interval check
641 }//p1 interval check
b966652f 642 }//end of 2nd particle loop
643 }//end of 1st particle loop
3feee083 644 //Printf("Number of analyzed events: %i",fAnalyzedEvents);
74b04dd4 645 //Printf("DeltaEta NN[0] = %.0f, PP[0] = %.0f, NP[0] = %.0f, PN[0] = %.0f",fNnn[kEta][0],fNpp[kEta][0],fNnp[kEta][0],fNpn[kEta][0]);
3feee083 646}
1d8792e7 647
1d8792e7 648
7f0257ea 649//____________________________________________________________________//
33a86282 650Double_t AliBalance::GetBalance(Int_t iAnalysisType, Int_t p2) {
6c178944 651 // Returns the value of the balance function in bin p2
33a86282 652 fB[iAnalysisType][p2] = 0.5*(((fNpn[iAnalysisType][p2] - 2.*fNnn[iAnalysisType][p2])/fNn[iAnalysisType]) + ((fNpn[iAnalysisType][p2] - 2.*fNpp[iAnalysisType][p2])/fNp[iAnalysisType]))/fP2Step[iAnalysisType];
7f0257ea 653
33a86282 654 return fB[iAnalysisType][p2];
6c178944 655}
656
7f0257ea 657//____________________________________________________________________//
33a86282 658Double_t AliBalance::GetError(Int_t iAnalysisType, Int_t p2) {
6c178944 659 // Returns the error on the BF value for bin p2
6432ac6a 660 // The errors for fNn and fNp are neglected here (0.1 % of total error)
33a86282 661 /*ferror[iAnalysisType][p2] = TMath::Sqrt(Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
662 + Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType]))
663 + Double_t(fNpn[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
664 + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
665 //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
666 //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3)
667 ) /fP2Step[iAnalysisType];*/
668
669 ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) +
670 Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) +
671 Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
6432ac6a 672
33a86282 673 return ferror[iAnalysisType][p2];
6c178944 674}
6432ac6a 675//____________________________________________________________________//
33a86282 676TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) {
6432ac6a 677
678 // Draws the BF
679 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
680 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
681 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
682 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
683
33a86282 684 if((fNp[iAnalysisType] == 0)||(fNn[iAnalysisType] == 0)) {
6432ac6a 685 cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
686 return NULL;
687 }
688
33a86282 689 for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
690 b[i] = GetBalance(iAnalysisType,i);
691 ber[i] = GetError(iAnalysisType,i);
692 x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
6432ac6a 693 xer[i] = 0.0;
694 }
695
33a86282 696 TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
6432ac6a 697 gr->GetXaxis()->SetTitleColor(1);
33a86282 698 if(iAnalysisType==0) {
6432ac6a 699 gr->SetTitle("Balance function B(#Delta y)");
700 gr->GetXaxis()->SetTitle("#Delta y");
701 gr->GetYaxis()->SetTitle("B(#Delta y)");
702 }
33a86282 703 if(iAnalysisType==1) {
6432ac6a 704 gr->SetTitle("Balance function B(#Delta #eta)");
705 gr->GetXaxis()->SetTitle("#Delta #eta");
706 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
707 }
33a86282 708 if(iAnalysisType==2) {
6432ac6a 709 gr->SetTitle("Balance function B(q_{long})");
710 gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
711 gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
712 }
33a86282 713 if(iAnalysisType==3) {
6432ac6a 714 gr->SetTitle("Balance function B(q_{out})");
715 gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
716 gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
717 }
33a86282 718 if(iAnalysisType==4) {
6432ac6a 719 gr->SetTitle("Balance function B(q_{side})");
720 gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
721 gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
722 }
33a86282 723 if(iAnalysisType==5) {
6432ac6a 724 gr->SetTitle("Balance function B(q_{inv})");
725 gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
726 gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
727 }
33a86282 728 if(iAnalysisType==6) {
6432ac6a 729 gr->SetTitle("Balance function B(#Delta #phi)");
730 gr->GetXaxis()->SetTitle("#Delta #phi");
731 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
732 }
33a86282 733
734 return gr;
6432ac6a 735}
6c178944 736
6432ac6a 737//____________________________________________________________________//
738void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
739 //Prints the calculated width of the BF and its error
740 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
741 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
742 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
743 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
744 Double_t deltaBalP2 = 0.0, integral = 0.0;
745 Double_t deltaErrorNew = 0.0;
746
747 cout<<"=================================================="<<endl;
748 for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
749 x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
7ef516b2 750 //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
6432ac6a 751 }
7ef516b2 752 //cout<<"=================================================="<<endl;
6432ac6a 753 for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
754 gSumXi += gHistBalance->GetBinCenter(i);
33a86282 755 gSumBi += gHistBalance->GetBinContent(i);
756 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
757 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
758 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
759 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
760 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
6432ac6a 761
33a86282 762 deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
763 integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
764 }
765 for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
766 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
767
768 Double_t integralError = TMath::Sqrt(deltaBalP2);
769
770 Double_t delta = gSumBiXi / gSumBi;
771 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
7ef516b2 772 cout<<"Analysis type: "<<gBFAnalysisType[iAnalysisType].Data()<<endl;
33a86282 773 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
774 cout<<"New error: "<<deltaErrorNew<<endl;
775 cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
776 cout<<"=================================================="<<endl;
6432ac6a 777}
778
779//____________________________________________________________________//
b7ce3d3f 780TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow) {
74b04dd4 781 //Returns the BF histogram, extracted from the 6 TH2D objects
6432ac6a 782 //(private members) of the AliBalance class.
b7ce3d3f 783 //
784 // Acceptance correction:
785 // - only for analysis type = kEta
786 // - only if etaWindow > 0 (default = -1.)
787 // - calculated as proposed by STAR
788 //
6432ac6a 789 TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
790 TString histName = "gHistBalanceFunctionHistogram";
791 histName += gAnalysisType[iAnalysisType];
792
74b04dd4 793 SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
794 fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
795 fHistPP[iAnalysisType]->GetNbinsY(),
796 fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),
797 fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
798
799 // determine the projection thresholds
800 Int_t binMinX, binMinY, binMinZ;
801 Int_t binMaxX, binMaxY, binMaxZ;
802
803 fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
804 fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);
6432ac6a 805
74b04dd4 806 TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
6432ac6a 807 switch(iAnalysisType) {
808 case kRapidity:
809 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
810 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
811 break;
812 case kEta:
813 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
814 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
815 break;
816 case kQlong:
817 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
818 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
819 break;
820 case kQout:
821 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
822 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
823 break;
824 case kQside:
825 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
826 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
827 break;
828 case kQinv:
829 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
830 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
831 break;
832 case kPhi:
833 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
834 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
835 break;
836 default:
837 break;
838 }
74b04dd4 839
840 TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
841 TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
842 TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
843 TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
844 TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
845 TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
33a86282 846
847 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
848 hTemp1->Sumw2();
849 hTemp2->Sumw2();
850 hTemp3->Sumw2();
851 hTemp4->Sumw2();
852 hTemp1->Add(hTemp3,-2.);
74b04dd4 853 hTemp1->Scale(1./hTemp5->GetEntries());
33a86282 854 hTemp2->Add(hTemp4,-2.);
74b04dd4 855 hTemp2->Scale(1./hTemp6->GetEntries());
33a86282 856 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
857 gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
858 }
b7ce3d3f 859
860 // do the acceptance correction (only for Eta and etaWindow > 0)
861 if(iAnalysisType == kEta && etaWindow > 0){
862 for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
863
864 Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
865 Double_t corrected = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
866 gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
867
868 }
869 }
6432ac6a 870
871 PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
898d5b28 872
6432ac6a 873 return gHistBalanceFunctionHistogram;
874}