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