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