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