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