1 /**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------
24 #include <Riostream.h>
30 #include <TLorentzVector.h>
31 #include <TObjArray.h>
32 #include <TGraphErrors.h>
35 #include "AliVParticle.h"
36 #include "AliMCParticle.h"
37 #include "AliESDtrack.h"
38 #include "AliAODTrack.h"
40 #include "AliBalance.h"
48 //____________________________________________________________________//
49 AliBalance::AliBalance() :
53 fConversionCut(kFALSE),
54 fAnalysisLevel("ESD"),
60 // Default constructor
62 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
64 fNumberOfBins[i] = 180;
72 fNumberOfBins[i] = 20;
78 fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
85 for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
104 fHistHBTbefore = NULL;
105 fHistHBTafter = NULL;
106 fHistConversionbefore = NULL;
107 fHistConversionafter = NULL;
112 //____________________________________________________________________//
113 AliBalance::AliBalance(const AliBalance& balance):
115 fShuffle(balance.fShuffle),
116 fHBTcut(balance.fHBTcut),
117 fConversionCut(balance.fConversionCut),
118 fAnalysisLevel(balance.fAnalysisLevel),
119 fAnalyzedEvents(balance.fAnalyzedEvents),
120 fCentralityId(balance.fCentralityId),
121 fCentStart(balance.fCentStart),
122 fCentStop(balance.fCentStop) {
124 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
125 fNn[i] = balance.fNn[i];
126 fNp[i] = balance.fNp[i];
128 fP1Start[i] = balance.fP1Start[i];
129 fP1Stop[i] = balance.fP1Stop[i];
130 fNumberOfBins[i] = balance.fNumberOfBins[i];
131 fP2Start[i] = balance.fP2Start[i];
132 fP2Stop[i] = balance.fP2Stop[i];
133 fP2Step[i] = balance.fP2Step[i];
134 fCentStart = balance.fCentStart;
135 fCentStop = balance.fCentStop;
137 fHistP[i] = balance.fHistP[i];
138 fHistN[i] = balance.fHistN[i];
139 fHistPN[i] = balance.fHistPN[i];
140 fHistNP[i] = balance.fHistNP[i];
141 fHistPP[i] = balance.fHistPP[i];
142 fHistNN[i] = balance.fHistNN[i];
144 for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
156 //____________________________________________________________________//
157 AliBalance::~AliBalance() {
160 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
172 //____________________________________________________________________//
173 void AliBalance::SetInterval(Int_t iAnalysisType,
174 Double_t p1Start, Double_t p1Stop,
175 Int_t ibins, Double_t p2Start, Double_t p2Stop) {
176 // Sets the analyzed interval.
177 // Set the same Information for all analyses
179 if(iAnalysisType == -1){
180 for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
181 fP1Start[i] = p1Start;
183 fNumberOfBins[i] = ibins;
184 fP2Start[i] = p2Start;
186 fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
189 // Set the Information for one analysis
190 else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
191 fP1Start[iAnalysisType] = p1Start;
192 fP1Stop[iAnalysisType] = p1Stop;
193 fNumberOfBins[iAnalysisType] = ibins;
194 fP2Start[iAnalysisType] = p2Start;
195 fP2Stop[iAnalysisType] = p2Stop;
196 fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
199 AliError("Wrong ANALYSIS number!");
203 //____________________________________________________________________//
204 void AliBalance::InitHistograms() {
205 //Initialize the histograms
207 // global switch disabling the reference
208 // (to avoid "Replacing existing TH1" if several wagons are created in train)
209 Bool_t oldStatus = TH1::AddDirectoryStatus();
210 TH1::AddDirectory(kFALSE);
213 for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
214 histName = "fHistP"; histName += kBFAnalysisType[iAnalysisType];
215 if(fShuffle) histName.Append("_shuffle");
216 if(fCentralityId) histName += fCentralityId.Data();
217 fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
219 histName = "fHistN"; histName += kBFAnalysisType[iAnalysisType];
220 if(fShuffle) histName.Append("_shuffle");
221 if(fCentralityId) histName += fCentralityId.Data();
222 fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
224 histName = "fHistPN"; histName += kBFAnalysisType[iAnalysisType];
225 if(fShuffle) histName.Append("_shuffle");
226 if(fCentralityId) histName += fCentralityId.Data();
227 fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
229 histName = "fHistNP"; histName += kBFAnalysisType[iAnalysisType];
230 if(fShuffle) histName.Append("_shuffle");
231 if(fCentralityId) histName += fCentralityId.Data();
232 fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
234 histName = "fHistPP"; histName += kBFAnalysisType[iAnalysisType];
235 if(fShuffle) histName.Append("_shuffle");
236 if(fCentralityId) histName += fCentralityId.Data();
237 fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
239 histName = "fHistNN"; histName += kBFAnalysisType[iAnalysisType];
240 if(fShuffle) histName.Append("_shuffle");
241 if(fCentralityId) histName += fCentralityId.Data();
242 fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
246 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200);
247 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200);
248 fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200);
249 fHistConversionafter = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200);
251 TH1::AddDirectory(oldStatus);
255 //____________________________________________________________________//
256 void AliBalance::PrintAnalysisSettings() {
257 //prints the analysis settings
259 Printf("======================================");
260 Printf("Analysis level: %s",fAnalysisLevel.Data());
261 Printf("======================================");
262 for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
263 Printf("Interval info for variable %d",ibin);
264 Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
265 Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
266 Printf("Number of bins: %d",fNumberOfBins[ibin]);
267 Printf("Step: %lf",fP2Step[ibin]);
270 Printf("======================================");
273 //____________________________________________________________________//
274 void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Float_t bSign) {
275 // Calculates the balance function
280 // Initialize histograms if not done yet
282 AliWarning("Histograms not yet initialized! --> Will be done now");
283 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
287 Int_t gNtrack = chargeVector[0]->size();
288 //Printf("(AliBalance) Number of tracks: %d",gNtrack);
290 for(i = 0; i < gNtrack;i++){
291 Short_t charge = chargeVector[0]->at(i);
292 Double_t rapidity = chargeVector[1]->at(i);
293 Double_t pseudorapidity = chargeVector[2]->at(i);
294 Double_t phi = chargeVector[3]->at(i);
296 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
297 for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
298 if(iAnalysisType == kEta) {
299 if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
301 fNp[iAnalysisType] += 1.;
302 fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
305 fNn[iAnalysisType] += 1.;
306 fHistN[iAnalysisType]->Fill(fCentrality,pseudorapidity);
309 }//analysis type: eta
310 else if(iAnalysisType == kPhi) {
311 if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
313 fNp[iAnalysisType] += 1.;
314 fHistP[iAnalysisType]->Fill(fCentrality,phi);
317 fNn[iAnalysisType] += 1.;
318 fHistN[iAnalysisType]->Fill(fCentrality,phi);
321 }//analysis type: phi
323 if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
325 fNp[iAnalysisType] += 1.;
326 fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
329 fNn[iAnalysisType] += 1.;
330 fHistN[iAnalysisType]->Fill(fCentrality,rapidity);
333 }//analysis type: y, qside, qout, qlong, qinv
334 }//analysis type loop
337 //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);
339 Double_t dy = 0., deta = 0.;
340 Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
344 Double_t eta1 = 0., rap1 = 0.;
345 Double_t px1 = 0., py1 = 0., pz1 = 0.;
347 Double_t energy1 = 0.;
351 Double_t eta2 = 0., rap2 = 0.;
352 Double_t px2 = 0., py2 = 0., pz2 = 0.;
354 Double_t energy2 = 0.;
356 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
357 for(i = 1; i < gNtrack; i++) {
359 charge1 = chargeVector[0]->at(i);
360 rap1 = chargeVector[1]->at(i);
361 eta1 = chargeVector[2]->at(i);
362 phi1 = chargeVector[3]->at(i);
363 px1 = chargeVector[4]->at(i);
364 py1 = chargeVector[5]->at(i);
365 pz1 = chargeVector[6]->at(i);
366 pt1 = chargeVector[7]->at(i);
367 energy1 = chargeVector[8]->at(i);
369 for(j = 0; j < i; j++) {
371 charge2 = chargeVector[0]->at(j);
372 rap2 = chargeVector[1]->at(j);
373 eta2 = chargeVector[2]->at(j);
374 phi2 = chargeVector[3]->at(j);
375 px2 = chargeVector[4]->at(j);
376 py2 = chargeVector[5]->at(j);
377 pz2 = chargeVector[6]->at(j);
378 pt2 = chargeVector[7]->at(j);
379 energy2 = chargeVector[8]->at(j);
381 // filling the arrays
384 dy = TMath::Abs(rap1 - rap2);
387 deta = TMath::Abs(eta1 - eta2);
390 Double_t eTot = energy1 + energy2;
391 Double_t pxTot = px1 + px2;
392 Double_t pyTot = py1 + py2;
393 Double_t pzTot = pz1 + pz2;
394 Double_t q0Tot = energy1 - energy2;
395 Double_t qxTot = px1 - px2;
396 Double_t qyTot = py1 - py2;
397 Double_t qzTot = pz1 - pz2;
399 Double_t eTot2 = eTot*eTot;
400 Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
401 Double_t pzTot2 = pzTot*pzTot;
403 Double_t q0Tot2 = q0Tot*q0Tot;
404 Double_t qTot2 = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;
406 Double_t snn = eTot2 - pTot2;
407 Double_t ptTot2 = pTot2 - pzTot2 ;
408 Double_t ptTot = TMath::Sqrt( ptTot2 );
410 qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
413 qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
416 qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
419 qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
422 dphi = TMath::Abs(phi1 - phi2);
423 if(dphi>180) dphi = 360 - dphi; //dphi should be between 0 and 180!
426 if(fHBTcut && charge1 * charge2 > 0){
427 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
430 // VERSION 2 (Taken from DPhiCorrelations)
431 // the variables & cuthave been developed by the HBT group
432 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
434 fHistHBTbefore->Fill(deta,dphi);
437 if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
441 Float_t phi1rad = phi1*TMath::DegToRad();
442 Float_t phi2rad = phi2*TMath::DegToRad();
444 // check first boundaries to see if is worth to loop and find the minimum
445 Float_t dphistar1 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 0.8, bSign);
446 Float_t dphistar2 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 2.5, bSign);
448 const Float_t kLimit = 0.02 * 3;
450 Float_t dphistarminabs = 1e5;
451 Float_t dphistarmin = 1e5;
453 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
455 for (Double_t rad=0.8; rad<2.51; rad+=0.01)
457 Float_t dphistar = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, rad, bSign);
458 Float_t dphistarabs = TMath::Abs(dphistar);
460 if (dphistarabs < dphistarminabs)
462 dphistarmin = dphistar;
463 dphistarminabs = dphistarabs;
467 if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
469 //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));
474 fHistHBTafter->Fill(deta,dphi);
479 if (charge1 * charge2 < 0)
482 fHistConversionbefore->Fill(deta,dphi);
484 Float_t m0 = 0.510e-3;
485 Float_t tantheta1 = 1e10;
488 Float_t phi1rad = phi1*TMath::DegToRad();
489 Float_t phi2rad = phi2*TMath::DegToRad();
491 if (eta1 < -1e-10 || eta1 > 1e-10)
492 tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
494 Float_t tantheta2 = 1e10;
495 if (eta2 < -1e-10 || eta2 > 1e-10)
496 tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
498 Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
499 Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
501 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
503 if (masssqu < 0.04*0.04){
504 //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));
507 fHistConversionafter->Fill(deta,dphi);
512 //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
513 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
516 if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
517 iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
518 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
520 if((charge1 > 0)&&(charge2 > 0)) {
521 fNpp[kRapidity][iBin] += 1.;
522 fHistPP[kRapidity]->Fill(fCentrality,dy);
524 else if((charge1 < 0)&&(charge2 < 0)) {
525 fNnn[kRapidity][iBin] += 1.;
526 fHistNN[kRapidity]->Fill(fCentrality,dy);
528 else if((charge1 > 0)&&(charge2 < 0)) {
529 fNpn[kRapidity][iBin] += 1.;
530 fHistPN[kRapidity]->Fill(fCentrality,dy);
532 else if((charge1 < 0)&&(charge2 > 0)) {
533 fNpn[kRapidity][iBin] += 1.;
534 fHistPN[kRapidity]->Fill(fCentrality,dy);
541 if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
542 if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
543 iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
544 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
545 if((charge1 > 0)&&(charge2 > 0)) {
546 fNpp[kEta][iBin] += 1.;
547 fHistPP[kEta]->Fill(fCentrality,deta);
549 if((charge1 < 0)&&(charge2 < 0)) {
550 fNnn[kEta][iBin] += 1.;
551 fHistNN[kEta]->Fill(fCentrality,deta);
553 if((charge1 > 0)&&(charge2 < 0)) {
554 fNpn[kEta][iBin] += 1.;
555 fHistPN[kEta]->Fill(fCentrality,deta);
557 if((charge1 < 0)&&(charge2 > 0)) {
558 fNpn[kEta][iBin] += 1.;
559 fHistPN[kEta]->Fill(fCentrality,deta);
565 // Qlong, out, side, inv
566 // Check the p1 intervall for rapidity here (like for single tracks above)
567 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
568 if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
569 iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
570 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
571 if((charge1 > 0)&&(charge2 > 0)) {
572 fNpp[kQlong][iBin] += 1.;
573 fHistPP[kQlong]->Fill(fCentrality,qLong);
575 if((charge1 < 0)&&(charge2 < 0)) {
576 fNnn[kQlong][iBin] += 1.;
577 fHistNN[kQlong]->Fill(fCentrality,qLong);
579 if((charge1 > 0)&&(charge2 < 0)) {
580 fNpn[kQlong][iBin] += 1.;
581 fHistPN[kQlong]->Fill(fCentrality,qLong);
583 if((charge1 < 0)&&(charge2 > 0)) {
584 fNpn[kQlong][iBin] += 1.;
585 fHistPN[kQlong]->Fill(fCentrality,qLong);
591 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
592 if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
593 iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
594 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
595 if((charge1 > 0)&&(charge2 > 0)) {
596 fNpp[kQout][iBin] += 1.;
597 fHistPP[kQout]->Fill(fCentrality,qOut);
599 if((charge1 < 0)&&(charge2 < 0)) {
600 fNnn[kQout][iBin] += 1.;
601 fHistNN[kQout]->Fill(fCentrality,qOut);
603 if((charge1 > 0)&&(charge2 < 0)) {
604 fNpn[kQout][iBin] += 1.;
605 fHistPN[kQout]->Fill(fCentrality,qOut);
607 if((charge1 < 0)&&(charge2 > 0)) {
608 fNpn[kQout][iBin] += 1.;
609 fHistPN[kQout]->Fill(fCentrality,qOut);
615 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
616 if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
617 iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
618 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
619 if((charge1 > 0)&&(charge2 > 0)) {
620 fNpp[kQside][iBin] += 1.;
621 fHistPP[kQside]->Fill(fCentrality,qSide);
623 if((charge1 < 0)&&(charge2 < 0)) {
624 fNnn[kQside][iBin] += 1.;
625 fHistNN[kQside]->Fill(fCentrality,qSide);
627 if((charge1 > 0)&&(charge2 < 0)) {
628 fNpn[kQside][iBin] += 1.;
629 fHistPN[kQside]->Fill(fCentrality,qSide);
631 if((charge1 < 0)&&(charge2 > 0)) {
632 fNpn[kQside][iBin] += 1.;
633 fHistPN[kQside]->Fill(fCentrality,qSide);
639 if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
640 if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
641 iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
642 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
643 if((charge1 > 0)&&(charge2 > 0)) {
644 fNpp[kQinv][iBin] += 1.;
645 fHistPP[kQinv]->Fill(fCentrality,qInv);
647 if((charge1 < 0)&&(charge2 < 0)) {
648 fNnn[kQinv][iBin] += 1.;
649 fHistNN[kQinv]->Fill(fCentrality,qInv);
651 if((charge1 > 0)&&(charge2 < 0)) {
652 fNpn[kQinv][iBin] += 1.;
653 fHistPN[kQinv]->Fill(fCentrality,qInv);
655 if((charge1 < 0)&&(charge2 > 0)) {
656 fNpn[kQinv][iBin] += 1.;
657 fHistPN[kQinv]->Fill(fCentrality,qInv);
664 if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
665 if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
666 iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
667 if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
668 if((charge1 > 0)&&(charge2 > 0)) {
669 fNpp[kPhi][iBin] += 1.;
670 fHistPP[kPhi]->Fill(fCentrality,dphi);
672 if((charge1 < 0)&&(charge2 < 0)) {
673 fNnn[kPhi][iBin] += 1.;
674 fHistNN[kPhi]->Fill(fCentrality,dphi);
676 if((charge1 > 0)&&(charge2 < 0)) {
677 fNpn[kPhi][iBin] += 1.;
678 fHistPN[kPhi]->Fill(fCentrality,dphi);
680 if((charge1 < 0)&&(charge2 > 0)) {
681 fNpn[kPhi][iBin] += 1.;
682 fHistPN[kPhi]->Fill(fCentrality,dphi);
687 }//end of 2nd particle loop
688 }//end of 1st particle loop
689 //Printf("Number of analyzed events: %i",fAnalyzedEvents);
690 //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]);
694 //____________________________________________________________________//
695 Double_t AliBalance::GetBalance(Int_t iAnalysisType, Int_t p2) {
696 // Returns the value of the balance function in bin p2
697 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];
699 return fB[iAnalysisType][p2];
702 //____________________________________________________________________//
703 Double_t AliBalance::GetError(Int_t iAnalysisType, Int_t p2) {
704 // Returns the error on the BF value for bin p2
705 // The errors for fNn and fNp are neglected here (0.1 % of total error)
706 /*ferror[iAnalysisType][p2] = TMath::Sqrt(Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
707 + Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType]))
708 + Double_t(fNpn[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
709 + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
710 //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
711 //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3)
712 ) /fP2Step[iAnalysisType];*/
714 ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) +
715 Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) +
716 Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
718 return ferror[iAnalysisType][p2];
720 //____________________________________________________________________//
721 TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) {
724 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
725 Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
726 Double_t b[MAXIMUM_NUMBER_OF_STEPS];
727 Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
729 if((fNp[iAnalysisType] == 0)||(fNn[iAnalysisType] == 0)) {
730 cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
734 for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
735 b[i] = GetBalance(iAnalysisType,i);
736 ber[i] = GetError(iAnalysisType,i);
737 x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
741 TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
742 gr->GetXaxis()->SetTitleColor(1);
743 if(iAnalysisType==0) {
744 gr->SetTitle("Balance function B(#Delta y)");
745 gr->GetXaxis()->SetTitle("#Delta y");
746 gr->GetYaxis()->SetTitle("B(#Delta y)");
748 if(iAnalysisType==1) {
749 gr->SetTitle("Balance function B(#Delta #eta)");
750 gr->GetXaxis()->SetTitle("#Delta #eta");
751 gr->GetYaxis()->SetTitle("B(#Delta #eta)");
753 if(iAnalysisType==2) {
754 gr->SetTitle("Balance function B(q_{long})");
755 gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
756 gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
758 if(iAnalysisType==3) {
759 gr->SetTitle("Balance function B(q_{out})");
760 gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
761 gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
763 if(iAnalysisType==4) {
764 gr->SetTitle("Balance function B(q_{side})");
765 gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
766 gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
768 if(iAnalysisType==5) {
769 gr->SetTitle("Balance function B(q_{inv})");
770 gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
771 gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
773 if(iAnalysisType==6) {
774 gr->SetTitle("Balance function B(#Delta #phi)");
775 gr->GetXaxis()->SetTitle("#Delta #phi");
776 gr->GetYaxis()->SetTitle("B(#Delta #phi)");
782 //____________________________________________________________________//
783 void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
784 //Prints the calculated width of the BF and its error
785 Double_t x[MAXIMUM_NUMBER_OF_STEPS];
786 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
787 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
788 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
789 Double_t deltaBalP2 = 0.0, integral = 0.0;
790 Double_t deltaErrorNew = 0.0;
792 cout<<"=================================================="<<endl;
793 for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) {
794 x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
795 //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
797 //cout<<"=================================================="<<endl;
798 for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
799 gSumXi += gHistBalance->GetBinCenter(i);
800 gSumBi += gHistBalance->GetBinContent(i);
801 gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
802 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
803 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
804 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
805 gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
807 deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
808 integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
810 for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
811 deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
813 Double_t integralError = TMath::Sqrt(deltaBalP2);
815 Double_t delta = gSumBiXi / gSumBi;
816 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
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;
824 //____________________________________________________________________//
825 TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly) {
826 //Returns the BF histogram, extracted from the 6 TH2D objects
827 //(private members) of the AliBalance class.
829 // Acceptance correction:
830 // - only for analysis type = kEta
831 // - only if etaWindow > 0 (default = -1.)
832 // - calculated as proposed by STAR
834 TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
835 TString histName = "gHistBalanceFunctionHistogram";
836 histName += gAnalysisType[iAnalysisType];
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());
844 // determine the projection thresholds
845 Int_t binMinX, binMinY, binMinZ;
846 Int_t binMaxX, binMaxY, binMaxZ;
848 fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
849 fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);
851 TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
852 switch(iAnalysisType) {
854 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
855 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
858 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
859 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
862 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
863 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
866 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
867 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
870 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
871 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
874 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
875 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
878 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
879 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
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));
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
895 TFile *fEfficiencyMatrix = NULL;
896 if(correctWithEfficiency){
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");
899 if(!fEfficiencyMatrix){
900 AliError("Efficiency histogram file not found");
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)
909 if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency){
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));
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));
922 // take the MC distributions
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));
928 if( !hEffP || !hEffN || !hEffPP || !hEffNN || !hEffPN){
929 AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
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))));
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))));
941 for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
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))));
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);
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)));
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)
977 if(iAnalysisType == kPhi && correctWithEfficiency){
979 TH1F* hEffPhiP = NULL;
980 TH1F* hEffPhiN = NULL;
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));
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));
990 // take the MC distributions
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));
996 if( !hEffPhiP || !hEffPhiN || !hEffPhiPP || !hEffPhiNN || !hEffPhiPN){
997 AliError("Efficiency (phi) histograms not found");
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))));
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))));
1009 for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
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))));
1025 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
1030 hTemp1->Add(hTemp3,-2.);
1031 hTemp1->Scale(1./hTemp5->Integral());
1032 hTemp2->Add(hTemp4,-2.);
1033 hTemp2->Scale(1./hTemp6->Integral());
1034 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1035 gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
1038 // do the acceptance correction (only for Eta and etaWindow > 0)
1039 if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency){
1040 for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
1042 Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
1043 Double_t corrected = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
1044 gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
1045 gHistBalanceFunctionHistogram->SetBinError(iBin+1,corrected/notCorrected*gHistBalanceFunctionHistogram->GetBinError(iBin+1));
1050 if(fEfficiencyMatrix) fEfficiencyMatrix->Close();
1052 PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);
1054 return gHistBalanceFunctionHistogram;