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