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