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 **************************************************************************/
14 /* $Id: AliBalancePsi.cxx 54125 2012-01-24 21:07:41Z miweber $ */
16 //-----------------------------------------------------------------
17 // Balance Function class
18 // This is the class to deal with the Balance Function wrt Psi analysis
19 // Origin: Panos Christakoglou, Nikhef, 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 "AliAnalysisTaskTriggeredBF.h"
42 #include "AliBalancePsi.h"
47 ClassImp(AliBalancePsi)
49 //____________________________________________________________________//
50 AliBalancePsi::AliBalancePsi() :
53 fAnalysisLevel("ESD"),
66 fHistConversionbefore(0),
67 fHistConversionafter(0),
69 fHistResonancesBefore(0),
70 fHistResonancesRho(0),
72 fHistResonancesLambda(0),
77 fResonancesCut(kFALSE),
80 fConversionCut(kFALSE),
81 fInvMassCutConversion(0.04),
84 fVertexBinning(kFALSE),
87 fEventClass("EventPlane"){
88 // Default constructor
91 //____________________________________________________________________//
92 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
93 TObject(balance), fShuffle(balance.fShuffle),
94 fAnalysisLevel(balance.fAnalysisLevel),
95 fAnalyzedEvents(balance.fAnalyzedEvents),
96 fCentralityId(balance.fCentralityId),
97 fCentStart(balance.fCentStart),
98 fCentStop(balance.fCentStop),
99 fHistP(balance.fHistP),
100 fHistN(balance.fHistN),
101 fHistPN(balance.fHistPN),
102 fHistNP(balance.fHistNP),
103 fHistPP(balance.fHistPP),
104 fHistNN(balance.fHistNN),
105 fHistHBTbefore(balance.fHistHBTbefore),
106 fHistHBTafter(balance.fHistHBTafter),
107 fHistConversionbefore(balance.fHistConversionbefore),
108 fHistConversionafter(balance.fHistConversionafter),
109 fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
110 fHistResonancesBefore(balance.fHistResonancesBefore),
111 fHistResonancesRho(balance.fHistResonancesRho),
112 fHistResonancesK0(balance.fHistResonancesK0),
113 fHistResonancesLambda(balance.fHistResonancesLambda),
114 fHistQbefore(balance.fHistQbefore),
115 fHistQafter(balance.fHistQafter),
116 fPsiInterval(balance.fPsiInterval),
117 fDeltaEtaMax(balance.fDeltaEtaMax),
118 fResonancesCut(balance.fResonancesCut),
119 fHBTCut(balance.fHBTCut),
120 fHBTCutValue(balance.fHBTCutValue),
121 fConversionCut(balance.fConversionCut),
122 fInvMassCutConversion(balance.fInvMassCutConversion),
123 fQCut(balance.fQCut),
124 fDeltaPtMin(balance.fDeltaPtMin),
125 fVertexBinning(balance.fVertexBinning),
126 fCustomBinning(balance.fCustomBinning),
127 fBinningString(balance.fBinningString),
128 fEventClass("EventPlane"){
132 //____________________________________________________________________//
133 AliBalancePsi::~AliBalancePsi() {
142 delete fHistHBTbefore;
143 delete fHistHBTafter;
144 delete fHistConversionbefore;
145 delete fHistConversionafter;
146 delete fHistPsiMinusPhi;
147 delete fHistResonancesBefore;
148 delete fHistResonancesRho;
149 delete fHistResonancesK0;
150 delete fHistResonancesLambda;
156 //____________________________________________________________________//
157 void AliBalancePsi::InitHistograms() {
158 // single particle histograms
160 // global switch disabling the reference
161 // (to avoid "Replacing existing TH1" if several wagons are created in train)
162 Bool_t oldStatus = TH1::AddDirectoryStatus();
163 TH1::AddDirectory(kFALSE);
165 Int_t anaSteps = 1; // analysis steps
166 Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
167 Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
168 TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
170 // two particle histograms
171 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
172 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
173 TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
177 // =========================================================
178 // The default string (from older versions of AliBalancePsi)
179 // =========================================================
180 TString defaultBinningStr;
181 defaultBinningStr = "multiplicity: 0,10,20,30,40,50,60,70,80,100,100000\n" // Multiplicity Bins
182 "centrality: 0.,5.,10.,20.,30.,40.,50.,60.,70.,80.\n" // Centrality Bins
183 "centralityVertex: 0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.\n" // Centrality Bins (Vertex Binning)
184 "eventPlane: -0.5,0.5,1.5,2.5,3.5\n" // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
185 "deltaEta: -1.6, -1.56, -1.52, -1.48, -1.44, -1.4, -1.36, -1.32, -1.28, -1.24, -1.2, -1.16, -1.12, -1.08, -1.04, -1, -0.96, -0.92, -0.88, -0.84, -0.8, -0.76, -0.72, -0.68, -0.64, -0.6, -0.56, -0.52, -0.48, -0.44, -0.4, -0.36, -0.32, -0.28, -0.24, -0.2, -0.16, -0.12, -0.08, -0.04, 0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6, 0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88, 0.92, 0.96, 1, 1.04, 1.08, 1.12, 1.16, 1.2, 1.24, 1.28, 1.32, 1.36, 1.4, 1.44, 1.48, 1.52, 1.56, 1.6\n" // Delta Eta Bins
186 "deltaEtaVertex: -1.6, -1.52, -1.44, -1.36, -1.28, -1.2, -1.12, -1.04, -0.96, -0.88, -0.8, -0.72, -0.64, -0.56, -0.48, -0.4, -0.32, -0.24, -0.16, -0.08, 0, 0.08, 0.16, 0.24, 0.32, 0.4, 0.48, 0.56, 0.64, 0.72, 0.8, 0.88, 0.96, 1.04, 1.12, 1.2, 1.28, 1.36, 1.44, 1.52, 1.6\n" // Delta Eta Bins (Vertex Binning)
187 "deltaPhi: -1.5708, -1.48353, -1.39626, -1.309, -1.22173, -1.13446, -1.0472, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.0872665, 0, 0.0872665, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.0472, 1.13446, 1.22173, 1.309, 1.39626, 1.48353, 1.5708, 1.65806, 1.74533, 1.8326, 1.91986, 2.00713, 2.0944, 2.18166, 2.26893, 2.35619, 2.44346, 2.53073, 2.61799, 2.70526, 2.79253, 2.87979, 2.96706, 3.05433, 3.14159, 3.22886, 3.31613, 3.40339, 3.49066, 3.57792, 3.66519, 3.75246, 3.83972, 3.92699, 4.01426, 4.10152, 4.18879, 4.27606, 4.36332, 4.45059, 4.53786, 4.62512, 4.71239\n" // Delta Phi Bins
188 "pT: 0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.\n" // pT Bins
189 "pTVertex: 0.2,1.0,2.0,3.0,4.0,8.0,15.0\n" // pT Bins (Vertex Binning)
190 "vertex: -10., 10.\n" // Vertex Bins
191 "vertexVertex: -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.\n" // Vertex Bins (Vertex Binning)
195 // =========================================================
196 // Customization (adopted from AliUEHistograms)
197 // =========================================================
199 TObjArray* lines = defaultBinningStr.Tokenize("\n");
200 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
202 TString line(lines->At(i)->GetName());
203 TString tag = line(0, line.Index(":")+1);
204 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
205 fBinningString += line + "\n";
207 AliInfo(Form("Using custom binning for %s", tag.Data()));
210 fBinningString += fCustomBinning;
212 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
215 // =========================================================
217 // =========================================================
219 //Depending on fEventClass Variable, do one thing or the other...
220 /**********************************************************
222 ======> Modification: Change Event Classification Scheme
224 ---> fEventClass == "EventPlane"
226 Default operation with Event Plane
228 ---> fEventClass == "Multiplicity"
230 Work with reference multiplicity (from GetReferenceMultiplicity, which one is decided in the configuration)
232 ---> fEventClass == "Centrality"
234 Work with Centrality Bins
236 ***********************************************************/
237 if(fEventClass == "Multiplicity"){
238 dBinsSingle[0] = GetBinning(fBinningString, "multiplicity", iBinSingle[0]);
239 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
240 axisTitleSingle[0] = "reference multiplicity";
241 axisTitlePair[0] = "reference multiplicity";
243 if(fEventClass == "Centrality"){
244 // fine binning in case of vertex Z binning
246 dBinsSingle[0] = GetBinning(fBinningString, "centralityVertex", iBinSingle[0]);
247 dBinsPair[0] = GetBinning(fBinningString, "centralityVertex", iBinPair[0]);
250 dBinsSingle[0] = GetBinning(fBinningString, "centrality", iBinSingle[0]);
251 dBinsPair[0] = GetBinning(fBinningString, "centrality", iBinPair[0]);
253 axisTitleSingle[0] = "Centrality percentile [%]";
254 axisTitlePair[0] = "Centrality percentile [%]";
256 if(fEventClass == "EventPlane"){
257 dBinsSingle[0] = GetBinning(fBinningString, "eventPlane", iBinSingle[0]);
258 dBinsPair[0] = GetBinning(fBinningString, "eventPlane", iBinPair[0]);
259 axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
260 axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
264 // Delta Eta and Delta Phi
265 // (coarse binning in case of vertex Z binning)
267 dBinsPair[1] = GetBinning(fBinningString, "deltaEtaVertex", iBinPair[1]);
270 dBinsPair[1] = GetBinning(fBinningString, "deltaEta", iBinPair[1]);
272 axisTitlePair[1] = "#Delta#eta";
274 dBinsPair[2] = GetBinning(fBinningString, "deltaPhi", iBinPair[2]);
275 axisTitlePair[2] = "#Delta#varphi (rad)";
278 // pT Trig and pT Assoc
279 // (coarse binning in case of vertex Z binning)
281 dBinsSingle[1] = GetBinning(fBinningString, "pTVertex", iBinSingle[1]);
282 dBinsPair[3] = GetBinning(fBinningString, "pTVertex", iBinPair[3]);
283 dBinsPair[4] = GetBinning(fBinningString, "pTVertex", iBinPair[4]);
286 dBinsSingle[1] = GetBinning(fBinningString, "pT", iBinSingle[1]);
287 dBinsPair[3] = GetBinning(fBinningString, "pT", iBinPair[3]);
288 dBinsPair[4] = GetBinning(fBinningString, "pT", iBinPair[4]);
291 axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
292 axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
293 axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
296 // vertex Z binning or not
298 dBinsSingle[2] = GetBinning(fBinningString, "vertexVertex", iBinSingle[2]);
299 dBinsPair[5] = GetBinning(fBinningString, "vertexVertex", iBinPair[5]);
302 dBinsSingle[2] = GetBinning(fBinningString, "vertex", iBinSingle[2]);
303 dBinsPair[5] = GetBinning(fBinningString, "vertex", iBinPair[5]);
306 axisTitleSingle[2] = "v_{Z} (cm)";
307 axisTitlePair[5] = "v_{Z} (cm)";
311 // =========================================================
312 // Create the Output objects (AliTHn)
313 // =========================================================
316 //+ triggered particles
318 if(fShuffle) histName.Append("_shuffle");
319 if(fCentralityId) histName += fCentralityId.Data();
320 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
321 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
322 fHistP->SetBinLimits(j, dBinsSingle[j]);
323 fHistP->SetVarTitle(j, axisTitleSingle[j]);
326 //- triggered particles
328 if(fShuffle) histName.Append("_shuffle");
329 if(fCentralityId) histName += fCentralityId.Data();
330 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
331 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
332 fHistN->SetBinLimits(j, dBinsSingle[j]);
333 fHistN->SetVarTitle(j, axisTitleSingle[j]);
337 histName = "fHistPN";
338 if(fShuffle) histName.Append("_shuffle");
339 if(fCentralityId) histName += fCentralityId.Data();
340 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
341 for (Int_t j=0; j<kTrackVariablesPair; j++) {
342 fHistPN->SetBinLimits(j, dBinsPair[j]);
343 fHistPN->SetVarTitle(j, axisTitlePair[j]);
347 histName = "fHistNP";
348 if(fShuffle) histName.Append("_shuffle");
349 if(fCentralityId) histName += fCentralityId.Data();
350 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
351 for (Int_t j=0; j<kTrackVariablesPair; j++) {
352 fHistNP->SetBinLimits(j, dBinsPair[j]);
353 fHistNP->SetVarTitle(j, axisTitlePair[j]);
357 histName = "fHistPP";
358 if(fShuffle) histName.Append("_shuffle");
359 if(fCentralityId) histName += fCentralityId.Data();
360 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
361 for (Int_t j=0; j<kTrackVariablesPair; j++) {
362 fHistPP->SetBinLimits(j, dBinsPair[j]);
363 fHistPP->SetVarTitle(j, axisTitlePair[j]);
367 histName = "fHistNN";
368 if(fShuffle) histName.Append("_shuffle");
369 if(fCentralityId) histName += fCentralityId.Data();
370 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
371 for (Int_t j=0; j<kTrackVariablesPair; j++) {
372 fHistNN->SetBinLimits(j, dBinsPair[j]);
373 fHistNN->SetVarTitle(j, axisTitlePair[j]);
376 AliInfo("Finished setting up the AliTHn");
379 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
380 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
381 fHistConversionbefore = new TH3D("fHistConversionbefore","before Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
382 fHistConversionafter = new TH3D("fHistConversionafter","after Conversion cut;#Delta#eta;#Delta#phi;M_{inv}^{2}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
383 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
384 fHistResonancesBefore = new TH3D("fHistResonancesBefore","before resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
385 fHistResonancesRho = new TH3D("fHistResonancesRho","after #rho resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
386 fHistResonancesK0 = new TH3D("fHistResonancesK0","after #rho, K0 resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
387 fHistResonancesLambda = new TH3D("fHistResonancesLambda","after #rho, K0, Lambda resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
388 fHistQbefore = new TH3D("fHistQbefore","before momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
389 fHistQafter = new TH3D("fHistQafter","after momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
391 TH1::AddDirectory(oldStatus);
395 //____________________________________________________________________//
396 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
397 TObjArray *particles,
398 TObjArray *particlesMixed,
400 Double_t kMultorCent,
402 // Calculates the balance function
405 // Initialize histograms if not done yet
407 AliWarning("Histograms not yet initialized! --> Will be done now");
408 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
412 Double_t trackVariablesSingle[kTrackVariablesSingle];
413 Double_t trackVariablesPair[kTrackVariablesPair];
416 AliWarning("particles TObjArray is NULL pointer --> return");
420 // define end of particle loops
421 Int_t iMax = particles->GetEntriesFast();
424 jMax = particlesMixed->GetEntriesFast();
426 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
427 TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
429 TArrayF secondEta(jMax);
430 TArrayF secondPhi(jMax);
431 TArrayF secondPt(jMax);
432 TArrayS secondCharge(jMax);
433 TArrayD secondCorrection(jMax);
435 for (Int_t i=0; i<jMax; i++){
436 secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
437 secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
438 secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
439 secondCharge[i] = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
440 secondCorrection[i] = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction(); //==========================correction
443 //TLorenzVector implementation for resonances
444 TLorentzVector vectorMother, vectorDaughter[2];
445 TParticle pPion, pProton, pRho0, pK0s, pLambda;
446 pPion.SetPdgCode(211); //pion
447 pRho0.SetPdgCode(113); //rho0
448 pK0s.SetPdgCode(310); //K0s
449 pProton.SetPdgCode(2212); //proton
450 pLambda.SetPdgCode(3122); //Lambda
451 Double_t gWidthForRho0 = 0.01;
452 Double_t gWidthForK0s = 0.01;
453 Double_t gWidthForLambda = 0.006;
454 Double_t nSigmaRejection = 3.0;
457 for (Int_t i = 0; i < iMax; i++) {
458 //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
459 AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
462 Float_t firstEta = firstParticle->Eta();
463 Float_t firstPhi = firstParticle->Phi();
464 Float_t firstPt = firstParticle->Pt();
465 Float_t firstCorrection = firstParticle->Correction();//==========================correction
467 // Event plane (determine psi bin)
468 Double_t gPsiMinusPhi = 0.;
469 Double_t gPsiMinusPhiBin = -10.;
470 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
472 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
473 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
474 gPsiMinusPhiBin = 0.0;
476 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
477 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
478 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
479 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
480 gPsiMinusPhiBin = 1.0;
482 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
483 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
484 gPsiMinusPhiBin = 2.0;
487 gPsiMinusPhiBin = 3.0;
489 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
491 Short_t charge1 = (Short_t) firstParticle->Charge();
493 trackVariablesSingle[0] = gPsiMinusPhiBin;
494 trackVariablesSingle[1] = firstPt;
495 if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
496 trackVariablesSingle[2] = vertexZ;
499 //fill single particle histograms
500 if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
501 else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
504 for(Int_t j = 0; j < jMax; j++) {
506 if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
508 // pT,Assoc < pT,Trig
509 if(firstPt < secondPt[j]) continue;
511 Short_t charge2 = secondCharge[j];
513 trackVariablesPair[0] = trackVariablesSingle[0];
514 trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
515 trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
516 //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
517 //trackVariablesPair[2] -= 360.;
518 //if (trackVariablesPair[2] < - 180.)
519 //trackVariablesPair[2] += 360.;
520 if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
521 trackVariablesPair[2] -= 2.*TMath::Pi();
522 if (trackVariablesPair[2] < - TMath::Pi())
523 trackVariablesPair[2] += 2.*TMath::Pi();
524 if (trackVariablesPair[2] < - TMath::Pi()/2.)
525 trackVariablesPair[2] += 2.*TMath::Pi();
527 trackVariablesPair[3] = firstPt; // pt trigger
528 trackVariablesPair[4] = secondPt[j]; // pt
529 trackVariablesPair[5] = vertexZ; // z of the primary vertex
531 //Exclude resonances for the calculation of pairs by looking
532 //at the invariant mass and not considering the pairs that
533 //fall within 3sigma from the mass peak of: rho0, K0s, Lambda
535 if (charge1 * charge2 < 0) {
538 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
539 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
540 vectorMother = vectorDaughter[0] + vectorDaughter[1];
541 fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
542 if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0)
544 fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
547 if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s)
549 fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
553 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
554 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass());
555 vectorMother = vectorDaughter[0] + vectorDaughter[1];
556 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
559 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass());
560 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
561 vectorMother = vectorDaughter[0] + vectorDaughter[1];
562 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
564 fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
570 if(fHBTCut){ // VERSION 3 (all pairs)
571 //if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
572 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
575 Double_t deta = firstEta - secondEta[j];
576 Double_t dphi = firstPhi - secondPhi[j];
577 // VERSION 2 (Taken from DPhiCorrelations)
578 // the variables & cuthave been developed by the HBT group
579 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
580 fHistHBTbefore->Fill(deta,dphi);
583 if (TMath::Abs(deta) < fHBTCutValue * 2.5 * 3) //fHBTCutValue = 0.02 [default for dphicorrelations]
586 //Float_t phi1rad = firstPhi*TMath::DegToRad();
587 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
588 Float_t phi1rad = firstPhi;
589 Float_t phi2rad = secondPhi[j];
591 // check first boundaries to see if is worth to loop and find the minimum
592 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
593 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
595 const Float_t kLimit = fHBTCutValue * 3;
597 Float_t dphistarminabs = 1e5;
598 Float_t dphistarmin = 1e5;
600 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
601 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
602 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
603 Float_t dphistarabs = TMath::Abs(dphistar);
605 if (dphistarabs < dphistarminabs) {
606 dphistarmin = dphistar;
607 dphistarminabs = dphistarabs;
611 if (dphistarminabs < fHBTCutValue && TMath::Abs(deta) < fHBTCutValue) {
612 //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));
617 fHistHBTafter->Fill(deta,dphi);
622 if (charge1 * charge2 < 0) {
623 Double_t deta = firstEta - secondEta[j];
624 Double_t dphi = firstPhi - secondPhi[j];
626 Float_t m0 = 0.510e-3;
627 Float_t tantheta1 = 1e10;
630 //Float_t phi1rad = firstPhi*TMath::DegToRad();
631 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
632 Float_t phi1rad = firstPhi;
633 Float_t phi2rad = secondPhi[j];
635 if (firstEta < -1e-10 || firstEta > 1e-10)
636 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
638 Float_t tantheta2 = 1e10;
639 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
640 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
642 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
643 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
645 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
647 fHistConversionbefore->Fill(deta,dphi,masssqu);
649 if (masssqu < fInvMassCutConversion*fInvMassCutConversion){
650 //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));
653 fHistConversionafter->Fill(deta,dphi,masssqu);
657 // momentum difference cut - suppress femtoscopic effects
660 //Double_t ptMin = 0.1; //const for the time being (should be changeable later on)
661 Double_t ptDifference = TMath::Abs( firstPt - secondPt[j]);
663 fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
664 if(ptDifference < fDeltaPtMin) continue;
665 fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
669 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
670 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
671 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
672 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
674 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
677 }//end of 2nd particle loop
678 }//end of 1st particle loop
681 //____________________________________________________________________//
682 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
688 Double_t ptTriggerMin,
689 Double_t ptTriggerMax,
690 Double_t ptAssociatedMin,
691 Double_t ptAssociatedMax) {
692 //Returns the BF histogram, extracted from the 6 AliTHn objects
693 //(private members) of the AliBalancePsi class.
694 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
695 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
698 if(psiMin > psiMax-0.00001){
699 AliError("psiMax <= psiMin");
702 if(vertexZMin > vertexZMax-0.00001){
703 AliError("vertexZMax <= vertexZMin");
706 if(ptTriggerMin > ptTriggerMax-0.00001){
707 AliError("ptTriggerMax <= ptTriggerMin");
710 if(ptAssociatedMin > ptAssociatedMax-0.00001){
711 AliError("ptAssociatedMax <= ptAssociatedMin");
716 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
717 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
718 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
719 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
720 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
721 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
725 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
726 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
727 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
728 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
729 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
730 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
734 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
735 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
736 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
737 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
738 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
739 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
740 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
744 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
745 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
746 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
747 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
748 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
751 //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
753 // Project into the wanted space (1st: analysis step, 2nd: axis)
754 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); //
755 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); //
756 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); //
757 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); //
758 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); //
759 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); //
761 TH1D *gHistBalanceFunctionHistogram = 0x0;
762 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
763 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
764 gHistBalanceFunctionHistogram->Reset();
766 switch(iVariablePair) {
768 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
769 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
772 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
773 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
783 hTemp1->Add(hTemp3,-1.);
784 hTemp1->Scale(1./hTemp5->Integral());
785 hTemp2->Add(hTemp4,-1.);
786 hTemp2->Scale(1./hTemp6->Integral());
787 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
788 gHistBalanceFunctionHistogram->Scale(0.5);
790 //normalize to bin width
791 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
794 return gHistBalanceFunctionHistogram;
797 //____________________________________________________________________//
798 TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
804 Double_t ptTriggerMin,
805 Double_t ptTriggerMax,
806 Double_t ptAssociatedMin,
807 Double_t ptAssociatedMax,
808 AliBalancePsi *bfMix) {
809 //Returns the BF histogram, extracted from the 6 AliTHn objects
810 //after dividing each correlation function by the Event Mixing one
811 //(private members) of the AliBalancePsi class.
812 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
813 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
816 if(psiMin > psiMax-0.00001){
817 AliError("psiMax <= psiMin");
820 if(vertexZMin > vertexZMax-0.00001){
821 AliError("vertexZMax <= vertexZMin");
824 if(ptTriggerMin > ptTriggerMax-0.00001){
825 AliError("ptTriggerMax <= ptTriggerMin");
828 if(ptAssociatedMin > ptAssociatedMax-0.00001){
829 AliError("ptAssociatedMax <= ptAssociatedMin");
833 // pt trigger (fixed for all vertices/psi/centralities)
834 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
835 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
836 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
837 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
838 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
839 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
840 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
843 // pt associated (fixed for all vertices/psi/centralities)
844 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
845 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
846 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
847 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
848 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
851 // ============================================================================================
852 // the same for event mixing
853 AliTHn *fHistPMix = bfMix->GetHistNp();
854 AliTHn *fHistNMix = bfMix->GetHistNn();
855 AliTHn *fHistPNMix = bfMix->GetHistNpn();
856 AliTHn *fHistNPMix = bfMix->GetHistNnp();
857 AliTHn *fHistPPMix = bfMix->GetHistNpp();
858 AliTHn *fHistNNMix = bfMix->GetHistNnn();
860 // pt trigger (fixed for all vertices/psi/centralities)
861 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
862 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
863 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
864 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
865 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
866 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
867 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
870 // pt associated (fixed for all vertices/psi/centralities)
871 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
872 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
873 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
874 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
875 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
878 // ============================================================================================
879 // ranges (in bins) for vertexZ and centrality (psi)
880 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
881 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
882 Int_t binVertexMin = 0;
883 Int_t binVertexMax = 0;
885 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
886 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
888 Double_t binPsiLowEdge = 0.;
889 Double_t binPsiUpEdge = 1000.;
890 Double_t binVertexLowEdge = 0.;
891 Double_t binVertexUpEdge = 1000.;
898 // loop over all bins
899 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
900 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
902 cout<<"In the balance function (1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
904 // determine the bin edges for this bin
905 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
906 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
908 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
909 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
913 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
914 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
915 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
916 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
917 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
918 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
922 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
923 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
924 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
925 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
926 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
927 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
930 // ============================================================================================
931 // the same for event mixing
934 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
935 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
936 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
937 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
938 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
939 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
943 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
944 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
945 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
946 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
947 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
948 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
951 // ============================================================================================
953 //Printf("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0));
955 // Project into the wanted space (1st: analysis step, 2nd: axis)
956 TH1D* hTempHelper1 = (TH1D*)fHistPN->Project(0,iVariablePair);
957 TH1D* hTempHelper2 = (TH1D*)fHistNP->Project(0,iVariablePair);
958 TH1D* hTempHelper3 = (TH1D*)fHistPP->Project(0,iVariablePair);
959 TH1D* hTempHelper4 = (TH1D*)fHistNN->Project(0,iVariablePair);
960 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
961 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
963 // ============================================================================================
964 // the same for event mixing
965 TH1D* hTempHelper1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
966 TH1D* hTempHelper2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
967 TH1D* hTempHelper3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
968 TH1D* hTempHelper4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
969 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
970 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
971 // ============================================================================================
973 hTempHelper1->Sumw2();
974 hTempHelper2->Sumw2();
975 hTempHelper3->Sumw2();
976 hTempHelper4->Sumw2();
979 hTempHelper1Mix->Sumw2();
980 hTempHelper2Mix->Sumw2();
981 hTempHelper3Mix->Sumw2();
982 hTempHelper4Mix->Sumw2();
986 // first put everything on positive x - axis (to be comparable with published data) --> ONLY IN THIS OPTION!
988 Double_t helperEndBin = 1.6;
989 if(iVariablePair==2) helperEndBin = TMath::Pi();
991 TH1D* hTempPos1 = new TH1D(Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
992 TH1D* hTempPos2 = new TH1D(Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
993 TH1D* hTempPos3 = new TH1D(Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
994 TH1D* hTempPos4 = new TH1D(Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
995 TH1D* hTempPos1Mix = new TH1D(Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
996 TH1D* hTempPos2Mix = new TH1D(Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
997 TH1D* hTempPos3Mix = new TH1D(Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
998 TH1D* hTempPos4Mix = new TH1D(Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1000 TH1D* hTemp1 = new TH1D(Form("hTemp1_%d_%d",iBinPsi,iBinVertex),Form("hTemp1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
1001 TH1D* hTemp2 = new TH1D(Form("hTemp2_%d_%d",iBinPsi,iBinVertex),Form("hTemp2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
1002 TH1D* hTemp3 = new TH1D(Form("hTemp3_%d_%d",iBinPsi,iBinVertex),Form("hTemp3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
1003 TH1D* hTemp4 = new TH1D(Form("hTemp4_%d_%d",iBinPsi,iBinVertex),Form("hTemp4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
1004 TH1D* hTemp1Mix = new TH1D(Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
1005 TH1D* hTemp2Mix = new TH1D(Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
1006 TH1D* hTemp3Mix = new TH1D(Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
1007 TH1D* hTemp4Mix = new TH1D(Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1014 hTempPos1Mix->Sumw2();
1015 hTempPos2Mix->Sumw2();
1016 hTempPos3Mix->Sumw2();
1017 hTempPos4Mix->Sumw2();
1029 for(Int_t i=0;i<hTempHelper1->GetNbinsX();i++){
1031 Double_t binCenter = hTempHelper1->GetXaxis()->GetBinCenter(i+1);
1033 if(iVariablePair==1){
1035 hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1036 hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1037 hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1038 hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1039 hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1040 hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1041 hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1042 hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1044 hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1045 hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1046 hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1047 hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1048 hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1049 hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1050 hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1051 hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1054 hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1055 hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1056 hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1057 hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1058 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1059 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1060 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1061 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1063 hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1064 hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1065 hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1066 hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1067 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1068 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1069 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1070 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1073 else if(iVariablePair==2){
1074 if(binCenter>0 && binCenter<TMath::Pi()){
1075 hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1076 hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1077 hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1078 hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1079 hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1080 hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1081 hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1082 hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1084 hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1085 hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1086 hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1087 hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1088 hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1089 hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1090 hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1091 hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1093 else if(binCenter<0){
1094 hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1095 hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1096 hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1097 hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1098 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1099 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1100 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1101 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1103 hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1104 hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1105 hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1106 hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1107 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1108 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1109 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1110 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1113 hTemp1->SetBinContent(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinContent(i+1));
1114 hTemp2->SetBinContent(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinContent(i+1));
1115 hTemp3->SetBinContent(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinContent(i+1));
1116 hTemp4->SetBinContent(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinContent(i+1));
1117 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1118 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1119 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1120 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1122 hTemp1->SetBinError(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinError(i+1));
1123 hTemp2->SetBinError(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinError(i+1));
1124 hTemp3->SetBinError(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinError(i+1));
1125 hTemp4->SetBinError(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinError(i+1));
1126 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinError(i+1));
1127 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinError(i+1));
1128 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinError(i+1));
1129 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinError(i+1));
1134 hTemp1->Add(hTempPos1);
1135 hTemp2->Add(hTempPos2);
1136 hTemp3->Add(hTempPos3);
1137 hTemp4->Add(hTempPos4);
1138 hTemp1Mix->Add(hTempPos1Mix);
1139 hTemp2Mix->Add(hTempPos2Mix);
1140 hTemp3Mix->Add(hTempPos3Mix);
1141 hTemp4Mix->Add(hTempPos4Mix);
1143 hTemp1->Scale(1./hTemp5->Integral());
1144 hTemp3->Scale(1./hTemp5->Integral());
1145 hTemp2->Scale(1./hTemp6->Integral());
1146 hTemp4->Scale(1./hTemp6->Integral());
1148 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1149 // does not work here, so normalize also to trigger particles
1150 // --> careful: gives different integrals then as with full 2D method
1151 hTemp1Mix->Scale(1./hTemp5Mix->Integral());
1152 hTemp3Mix->Scale(1./hTemp5Mix->Integral());
1153 hTemp2Mix->Scale(1./hTemp6Mix->Integral());
1154 hTemp4Mix->Scale(1./hTemp6Mix->Integral());
1156 hTemp1->Divide(hTemp1Mix);
1157 hTemp2->Divide(hTemp2Mix);
1158 hTemp3->Divide(hTemp3Mix);
1159 hTemp4->Divide(hTemp4Mix);
1161 // for the first: clone
1162 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1163 h1 = (TH1D*)hTemp1->Clone();
1164 h2 = (TH1D*)hTemp2->Clone();
1165 h3 = (TH1D*)hTemp3->Clone();
1166 h4 = (TH1D*)hTemp4->Clone();
1168 else{ // otherwise: add for averaging
1188 delete hTempPos1Mix;
1189 delete hTempPos2Mix;
1190 delete hTempPos3Mix;
1191 delete hTempPos4Mix;
1195 TH1D *gHistBalanceFunctionHistogram = 0x0;
1196 if((h1)&&(h2)&&(h3)&&(h4)) {
1198 // average over number of bins nbinsVertex * nbinsPsi
1199 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1200 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1201 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1202 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1204 gHistBalanceFunctionHistogram = (TH1D*)h1->Clone();
1205 gHistBalanceFunctionHistogram->Reset();
1207 switch(iVariablePair) {
1209 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1210 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1213 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1214 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1223 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1224 gHistBalanceFunctionHistogram->Scale(0.5);
1226 //normalize to bin width
1227 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1230 return gHistBalanceFunctionHistogram;
1233 //____________________________________________________________________//
1234 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
1236 Double_t vertexZMin,
1237 Double_t vertexZMax,
1238 Double_t ptTriggerMin,
1239 Double_t ptTriggerMax,
1240 Double_t ptAssociatedMin,
1241 Double_t ptAssociatedMax) {
1242 //Returns the BF histogram in Delta eta vs Delta phi,
1243 //extracted from the 6 AliTHn objects
1244 //(private members) of the AliBalancePsi class.
1245 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1246 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1249 if(psiMin > psiMax-0.00001){
1250 AliError("psiMax <= psiMin");
1253 if(vertexZMin > vertexZMax-0.00001){
1254 AliError("vertexZMax <= vertexZMin");
1257 if(ptTriggerMin > ptTriggerMax-0.00001){
1258 AliError("ptTriggerMax <= ptTriggerMin");
1261 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1262 AliError("ptAssociatedMax <= ptAssociatedMin");
1266 TString histName = "gHistBalanceFunctionHistogram2D";
1269 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1270 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1271 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1272 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1273 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1274 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1278 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1279 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1280 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1281 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1282 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1283 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1287 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1288 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1289 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1290 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1291 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1292 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1293 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1297 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1298 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1299 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1300 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1301 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1304 //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
1306 // Project into the wanted space (1st: analysis step, 2nd: axis)
1307 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1308 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1309 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1310 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1311 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1312 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1314 TH2D *gHistBalanceFunctionHistogram = 0x0;
1315 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1316 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1317 gHistBalanceFunctionHistogram->Reset();
1318 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1319 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1320 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1326 hTemp1->Add(hTemp3,-1.);
1327 hTemp1->Scale(1./hTemp5->Integral());
1328 hTemp2->Add(hTemp4,-1.);
1329 hTemp2->Scale(1./hTemp6->Integral());
1330 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1331 gHistBalanceFunctionHistogram->Scale(0.5);
1333 //normalize to bin width
1334 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1337 return gHistBalanceFunctionHistogram;
1340 //____________________________________________________________________//
1341 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1343 Double_t vertexZMin,
1344 Double_t vertexZMax,
1345 Double_t ptTriggerMin,
1346 Double_t ptTriggerMax,
1347 Double_t ptAssociatedMin,
1348 Double_t ptAssociatedMax,
1349 AliBalancePsi *bfMix) {
1350 //Returns the BF histogram in Delta eta vs Delta phi,
1351 //after dividing each correlation function by the Event Mixing one
1352 //extracted from the 6 AliTHn objects
1353 //(private members) of the AliBalancePsi class.
1354 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1355 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1356 TString histName = "gHistBalanceFunctionHistogram2D";
1359 AliError("balance function object for event mixing not available");
1364 if(psiMin > psiMax-0.00001){
1365 AliError("psiMax <= psiMin");
1368 if(vertexZMin > vertexZMax-0.00001){
1369 AliError("vertexZMax <= vertexZMin");
1372 if(ptTriggerMin > ptTriggerMax-0.00001){
1373 AliError("ptTriggerMax <= ptTriggerMin");
1376 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1377 AliError("ptAssociatedMax <= ptAssociatedMin");
1381 // pt trigger (fixed for all vertices/psi/centralities)
1382 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1383 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1384 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1385 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1386 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1387 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1388 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1391 // pt associated (fixed for all vertices/psi/centralities)
1392 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1393 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1394 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1395 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1396 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1399 // ============================================================================================
1400 // the same for event mixing
1401 AliTHn *fHistPMix = bfMix->GetHistNp();
1402 AliTHn *fHistNMix = bfMix->GetHistNn();
1403 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1404 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1405 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1406 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1408 // pt trigger (fixed for all vertices/psi/centralities)
1409 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1410 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1411 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1412 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1413 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1414 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1415 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1418 // pt associated (fixed for all vertices/psi/centralities)
1419 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1420 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1421 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1422 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1423 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1426 // ============================================================================================
1427 // ranges (in bins) for vertexZ and centrality (psi)
1428 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1429 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1430 Int_t binVertexMin = 0;
1431 Int_t binVertexMax = 0;
1433 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1434 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1436 Double_t binPsiLowEdge = 0.;
1437 Double_t binPsiUpEdge = 0.;
1438 Double_t binVertexLowEdge = 0.;
1439 Double_t binVertexUpEdge = 0.;
1446 // loop over all bins
1447 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1448 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1450 cout<<"In the balance function (2D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1452 // determine the bin edges for this bin
1453 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1454 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1456 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1457 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1461 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1462 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1463 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1464 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1465 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1466 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1470 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1471 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1472 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1473 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1474 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1475 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1480 // ============================================================================================
1481 // the same for event mixing
1484 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1485 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1486 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1487 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1488 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1489 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1493 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1494 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1495 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1496 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1497 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1498 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1501 // ============================================================================================
1504 //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
1506 // Project into the wanted space (1st: analysis step, 2nd: axis)
1507 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1508 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1509 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1510 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1511 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1512 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1514 // ============================================================================================
1515 // the same for event mixing
1516 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1517 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1518 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1519 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1520 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1521 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1522 // ============================================================================================
1533 hTemp1->Scale(1./hTemp5->Integral());
1534 hTemp3->Scale(1./hTemp5->Integral());
1535 hTemp2->Scale(1./hTemp6->Integral());
1536 hTemp4->Scale(1./hTemp6->Integral());
1538 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1539 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1540 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1541 hTemp1Mix->Scale(1./mixedNorm1);
1542 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1543 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1544 hTemp2Mix->Scale(1./mixedNorm2);
1545 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1546 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1547 hTemp3Mix->Scale(1./mixedNorm3);
1548 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1549 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1550 hTemp4Mix->Scale(1./mixedNorm4);
1552 hTemp1->Divide(hTemp1Mix);
1553 hTemp2->Divide(hTemp2Mix);
1554 hTemp3->Divide(hTemp3Mix);
1555 hTemp4->Divide(hTemp4Mix);
1557 // for the first: clone
1558 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1559 h1 = (TH2D*)hTemp1->Clone();
1560 h2 = (TH2D*)hTemp2->Clone();
1561 h3 = (TH2D*)hTemp3->Clone();
1562 h4 = (TH2D*)hTemp4->Clone();
1564 else{ // otherwise: add for averaging
1573 TH2D *gHistBalanceFunctionHistogram = 0x0;
1574 if((h1)&&(h2)&&(h3)&&(h4)) {
1576 // average over number of bins nbinsVertex * nbinsPsi
1577 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1578 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1579 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1580 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1582 gHistBalanceFunctionHistogram = (TH2D*)h1->Clone();
1583 gHistBalanceFunctionHistogram->Reset();
1584 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1585 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1586 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1591 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1592 gHistBalanceFunctionHistogram->Scale(0.5);
1594 //normalize to bin width
1595 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1598 return gHistBalanceFunctionHistogram;
1601 //____________________________________________________________________//
1602 TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1605 Double_t vertexZMin,
1606 Double_t vertexZMax,
1607 Double_t ptTriggerMin,
1608 Double_t ptTriggerMax,
1609 Double_t ptAssociatedMin,
1610 Double_t ptAssociatedMax,
1611 AliBalancePsi *bfMix) {
1612 //Returns the BF histogram in Delta eta OR Delta phi,
1613 //after dividing each correlation function by the Event Mixing one
1614 // (But the division is done here in 2D, this was basically done to check the Integral)
1615 //extracted from the 6 AliTHn objects
1616 //(private members) of the AliBalancePsi class.
1617 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1618 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1619 TString histName = "gHistBalanceFunctionHistogram1D";
1622 AliError("balance function object for event mixing not available");
1627 if(psiMin > psiMax-0.00001){
1628 AliError("psiMax <= psiMin");
1631 if(vertexZMin > vertexZMax-0.00001){
1632 AliError("vertexZMax <= vertexZMin");
1635 if(ptTriggerMin > ptTriggerMax-0.00001){
1636 AliError("ptTriggerMax <= ptTriggerMin");
1639 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1640 AliError("ptAssociatedMax <= ptAssociatedMin");
1644 // pt trigger (fixed for all vertices/psi/centralities)
1645 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1646 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1647 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1648 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1649 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1650 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1651 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1654 // pt associated (fixed for all vertices/psi/centralities)
1655 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1656 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1657 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1658 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1659 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1662 // ============================================================================================
1663 // the same for event mixing
1664 AliTHn *fHistPMix = bfMix->GetHistNp();
1665 AliTHn *fHistNMix = bfMix->GetHistNn();
1666 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1667 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1668 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1669 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1671 // pt trigger (fixed for all vertices/psi/centralities)
1672 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1673 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1674 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1675 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1676 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1677 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1678 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1681 // pt associated (fixed for all vertices/psi/centralities)
1682 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1683 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1684 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1685 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1686 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1689 // ============================================================================================
1690 // ranges (in bins) for vertexZ and centrality (psi)
1691 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1692 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1693 Int_t binVertexMin = 0;
1694 Int_t binVertexMax = 0;
1696 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1697 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1699 Double_t binPsiLowEdge = 0.;
1700 Double_t binPsiUpEdge = 0.;
1701 Double_t binVertexLowEdge = 0.;
1702 Double_t binVertexUpEdge = 0.;
1709 // loop over all bins
1710 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1711 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1713 cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1715 // determine the bin edges for this bin
1716 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1717 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1719 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1720 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1724 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1725 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1726 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1727 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1728 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1729 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1733 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1734 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1735 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1736 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1737 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1738 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1741 // ============================================================================================
1742 // the same for event mixing
1745 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1746 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1747 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1748 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1749 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1750 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1754 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1755 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1756 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1757 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1758 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1759 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1761 // ============================================================================================
1763 //AliInfo(Form("P:%lf - N:%lf - PN:%lf - NP:%lf - PP:%lf - NN:%lf",fHistP->GetEntries(0),fHistN->GetEntries(0),fHistPN->GetEntries(0),fHistNP->GetEntries(0),fHistPP->GetEntries(0),fHistNN->GetEntries(0)));
1765 // Project into the wanted space (1st: analysis step, 2nd: axis)
1766 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1767 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1768 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1769 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1770 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1771 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1773 // ============================================================================================
1774 // the same for event mixing
1775 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1776 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1777 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1778 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1779 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1780 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1781 // ============================================================================================
1792 hTemp1->Scale(1./hTemp5->Integral());
1793 hTemp3->Scale(1./hTemp5->Integral());
1794 hTemp2->Scale(1./hTemp6->Integral());
1795 hTemp4->Scale(1./hTemp6->Integral());
1797 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1798 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1799 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1800 hTemp1Mix->Scale(1./mixedNorm1);
1801 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1802 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1803 hTemp2Mix->Scale(1./mixedNorm2);
1804 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1805 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1806 hTemp3Mix->Scale(1./mixedNorm3);
1807 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1808 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1809 hTemp4Mix->Scale(1./mixedNorm4);
1811 hTemp1->Divide(hTemp1Mix);
1812 hTemp2->Divide(hTemp2Mix);
1813 hTemp3->Divide(hTemp3Mix);
1814 hTemp4->Divide(hTemp4Mix);
1816 // for the first: clone
1817 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1818 h1 = (TH2D*)hTemp1->Clone();
1819 h2 = (TH2D*)hTemp2->Clone();
1820 h3 = (TH2D*)hTemp3->Clone();
1821 h4 = (TH2D*)hTemp4->Clone();
1823 else{ // otherwise: add for averaging
1833 TH1D *gHistBalanceFunctionHistogram = 0x0;
1834 if((h1)&&(h2)&&(h3)&&(h4)) {
1836 // average over number of bins nbinsVertex * nbinsPsi
1837 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1838 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1839 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1840 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1842 // now only project on one axis
1843 TH1D *h1DTemp1 = NULL;
1844 TH1D *h1DTemp2 = NULL;
1845 TH1D *h1DTemp3 = NULL;
1846 TH1D *h1DTemp4 = NULL;
1848 h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName()));
1849 h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName()));
1850 h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName()));
1851 h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName()));
1854 h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName()));
1855 h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName()));
1856 h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName()));
1857 h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName()));
1860 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1861 gHistBalanceFunctionHistogram->Reset();
1863 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1864 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1867 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1868 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1871 h1DTemp1->Add(h1DTemp3,-1.);
1872 h1DTemp2->Add(h1DTemp4,-1.);
1874 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1875 gHistBalanceFunctionHistogram->Scale(0.5);
1877 //normalize to bin width
1878 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1881 return gHistBalanceFunctionHistogram;
1885 //____________________________________________________________________//
1886 TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1889 Double_t vertexZMin,
1890 Double_t vertexZMax,
1891 Double_t ptTriggerMin,
1892 Double_t ptTriggerMax,
1893 Double_t ptAssociatedMin,
1894 Double_t ptAssociatedMax,
1895 AliBalancePsi *bMixed,
1896 Bool_t normToTrig) {
1898 // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1899 // does the division by event mixing inside,
1900 // and averaging over several vertexZ and centrality bins
1903 if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1904 AliError("Only types allowed: PN,NP,PP,NN,ALL");
1908 AliError("No Event Mixing AliTHn");
1914 TH2D *fMixed = NULL;
1916 // ranges (in bins) for vertexZ and centrality (psi)
1917 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1918 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1919 Int_t binVertexMin = 0;
1920 Int_t binVertexMax = 0;
1922 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1923 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1925 Double_t binPsiLowEdge = 0.;
1926 Double_t binPsiUpEdge = 1000.;
1927 Double_t binVertexLowEdge = 0.;
1928 Double_t binVertexUpEdge = 1000.;
1930 // loop over all bins
1931 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1932 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1934 cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1936 // determine the bin edges for this bin
1937 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1938 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1940 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1941 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1944 // get the 2D histograms for the correct type
1946 fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1947 fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1949 else if(type=="NP"){
1950 fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1951 fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1953 else if(type=="PP"){
1954 fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1955 fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1957 else if(type=="NN"){
1958 fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1959 fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1961 else if(type=="ALL"){
1962 fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1963 fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1966 if(fMixed && normToTrig && fMixed->Integral()>0){
1968 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1969 // do it only on away-side (due to two-track cuts)
1970 Double_t mixedNorm = fMixed->Integral(fMixed->GetXaxis()->FindBin(0-10e-5),fMixed->GetXaxis()->FindBin(0+10e-5),fMixed->GetNbinsY()/2+1,fMixed->GetNbinsY());
1971 mixedNorm /= 0.5 * fMixed->GetNbinsY() *(fMixed->GetXaxis()->FindBin(0.01) - fMixed->GetXaxis()->FindBin(-0.01) + 1);
1973 // finite bin correction
1974 Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
1975 Double_t maxEta = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
1977 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1978 //Printf("Finite bin correction: %f", finiteBinCorrection);
1979 mixedNorm /= finiteBinCorrection;
1981 fMixed->Scale(1./mixedNorm);
1984 if(fSame && fMixed){
1985 // then get the correlation function (divide fSame/fmixed)
1986 fSame->Divide(fMixed);
1989 // average over number of triggers in each sub-bin
1990 Double_t NTrigSubBin = 0;
1991 if(type=="PN" || type=="PP")
1992 NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
1993 else if(type=="NP" || type=="NN")
1994 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
1995 fSame->Scale(NTrigSubBin);
1997 // for the first: clone
1998 if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
1999 gHist = (TH2D*)fSame->Clone();
2001 else{ // otherwise: add for averaging
2011 // average over number of bins nbinsVertex * nbinsPsi
2012 // gHist->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
2015 // average over number of triggers in each sub-bin
2016 // first set to full range and then obtain number of all triggers
2017 Double_t NTrigAll = 0;
2018 if(type=="PN" || type=="PP"){
2019 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2020 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2021 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2022 NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
2024 else if(type=="NP" || type=="NN"){
2025 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2026 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2027 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2028 NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
2030 gHist->Scale(1./NTrigAll);
2038 //____________________________________________________________________//
2039 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
2041 Double_t vertexZMin,
2042 Double_t vertexZMax,
2043 Double_t ptTriggerMin,
2044 Double_t ptTriggerMax,
2045 Double_t ptAssociatedMin,
2046 Double_t ptAssociatedMax) {
2047 //Returns the 2D correlation function for (+-) pairs
2050 if(psiMin > psiMax-0.00001){
2051 AliError("psiMax <= psiMin");
2054 if(vertexZMin > vertexZMax-0.00001){
2055 AliError("vertexZMax <= vertexZMin");
2058 if(ptTriggerMin > ptTriggerMax-0.00001){
2059 AliError("ptTriggerMax <= ptTriggerMin");
2062 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2063 AliError("ptAssociatedMax <= ptAssociatedMin");
2068 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2069 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2073 //Printf("Cutting on Vz...");
2074 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2075 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2079 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2080 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2081 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2085 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2086 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2088 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2089 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2091 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
2092 //TCanvas *c1 = new TCanvas("c1","");
2095 //AliError("Projection of fHistP = NULL");
2099 //gHistTest->DrawCopy("colz");
2102 //0:step, 1: Delta eta, 2: Delta phi
2103 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2105 AliError("Projection of fHistPN = NULL");
2109 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
2110 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
2111 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
2112 //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
2113 //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
2115 //TCanvas *c2 = new TCanvas("c2","");
2117 //fHistPN->Project(0,1,2)->DrawCopy("colz");
2119 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2120 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2122 //normalize to bin width
2123 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2128 //____________________________________________________________________//
2129 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
2131 Double_t vertexZMin,
2132 Double_t vertexZMax,
2133 Double_t ptTriggerMin,
2134 Double_t ptTriggerMax,
2135 Double_t ptAssociatedMin,
2136 Double_t ptAssociatedMax) {
2137 //Returns the 2D correlation function for (+-) pairs
2140 if(psiMin > psiMax-0.00001){
2141 AliError("psiMax <= psiMin");
2144 if(vertexZMin > vertexZMax-0.00001){
2145 AliError("vertexZMax <= vertexZMin");
2148 if(ptTriggerMin > ptTriggerMax-0.00001){
2149 AliError("ptTriggerMax <= ptTriggerMin");
2152 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2153 AliError("ptAssociatedMax <= ptAssociatedMin");
2158 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2159 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2163 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2164 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2168 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2169 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2170 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2174 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2175 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2177 //0:step, 1: Delta eta, 2: Delta phi
2178 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2180 AliError("Projection of fHistPN = NULL");
2184 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2185 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
2186 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2187 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2189 //normalize to bin width
2190 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2195 //____________________________________________________________________//
2196 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
2198 Double_t vertexZMin,
2199 Double_t vertexZMax,
2200 Double_t ptTriggerMin,
2201 Double_t ptTriggerMax,
2202 Double_t ptAssociatedMin,
2203 Double_t ptAssociatedMax) {
2204 //Returns the 2D correlation function for (+-) pairs
2207 if(psiMin > psiMax-0.00001){
2208 AliError("psiMax <= psiMin");
2211 if(vertexZMin > vertexZMax-0.00001){
2212 AliError("vertexZMax <= vertexZMin");
2215 if(ptTriggerMin > ptTriggerMax-0.00001){
2216 AliError("ptTriggerMax <= ptTriggerMin");
2219 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2220 AliError("ptAssociatedMax <= ptAssociatedMin");
2225 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2226 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2230 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2231 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2235 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2236 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2237 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2241 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2242 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2244 //0:step, 1: Delta eta, 2: Delta phi
2245 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2247 AliError("Projection of fHistPN = NULL");
2251 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2252 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
2253 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2254 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2256 //normalize to bin width
2257 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2262 //____________________________________________________________________//
2263 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
2265 Double_t vertexZMin,
2266 Double_t vertexZMax,
2267 Double_t ptTriggerMin,
2268 Double_t ptTriggerMax,
2269 Double_t ptAssociatedMin,
2270 Double_t ptAssociatedMax) {
2271 //Returns the 2D correlation function for (+-) pairs
2274 if(psiMin > psiMax-0.00001){
2275 AliError("psiMax <= psiMin");
2278 if(vertexZMin > vertexZMax-0.00001){
2279 AliError("vertexZMax <= vertexZMin");
2282 if(ptTriggerMin > ptTriggerMax-0.00001){
2283 AliError("ptTriggerMax <= ptTriggerMin");
2286 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2287 AliError("ptAssociatedMax <= ptAssociatedMin");
2292 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2293 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2297 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2298 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2302 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2303 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2304 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2308 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2309 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2311 //0:step, 1: Delta eta, 2: Delta phi
2312 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2314 AliError("Projection of fHistPN = NULL");
2318 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2319 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
2320 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2321 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2323 //normalize to bin width
2324 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2329 //____________________________________________________________________//
2330 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2332 Double_t vertexZMin,
2333 Double_t vertexZMax,
2334 Double_t ptTriggerMin,
2335 Double_t ptTriggerMax,
2336 Double_t ptAssociatedMin,
2337 Double_t ptAssociatedMax) {
2338 //Returns the 2D correlation function for the sum of all charge combination pairs
2341 if(psiMin > psiMax-0.00001){
2342 AliError("psiMax <= psiMin");
2345 if(vertexZMin > vertexZMax-0.00001){
2346 AliError("vertexZMax <= vertexZMin");
2349 if(ptTriggerMin > ptTriggerMax-0.00001){
2350 AliError("ptTriggerMax <= ptTriggerMin");
2353 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2354 AliError("ptAssociatedMax <= ptAssociatedMin");
2359 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2360 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2361 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2362 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2363 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2364 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2368 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2369 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2370 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2371 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2372 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2373 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2377 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2378 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2379 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2380 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2381 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2382 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2383 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2387 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2388 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2389 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2390 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2391 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2393 //0:step, 1: Delta eta, 2: Delta phi
2394 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2396 AliError("Projection of fHistNN = NULL");
2399 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2401 AliError("Projection of fHistPP = NULL");
2404 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2406 AliError("Projection of fHistNP = NULL");
2409 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2411 AliError("Projection of fHistPN = NULL");
2415 // sum all 2 particle histograms
2416 gHistNN->Add(gHistPP);
2417 gHistNN->Add(gHistNP);
2418 gHistNN->Add(gHistPN);
2420 // divide by sum of + and - triggers
2421 if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
2422 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistN->Project(0,1)->Integral()));
2424 //normalize to bin width
2425 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
2430 //____________________________________________________________________//
2432 Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
2433 Double_t &mean, Double_t &meanError,
2434 Double_t &sigma, Double_t &sigmaError,
2435 Double_t &skewness, Double_t &skewnessError,
2436 Double_t &kurtosis, Double_t &kurtosisError) {
2438 // helper method to calculate the moments and errors of a TH1D anlytically
2439 // fVariable = 1 for Delta eta (moments in full range)
2440 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
2443 Bool_t success = kFALSE;
2449 skewnessError = -1.;
2451 kurtosisError = -1.;
2455 // ----------------------------------------------------------------------
2456 // basic parameters of histogram
2458 Int_t fNumberOfBins = gHist->GetNbinsX();
2459 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
2462 // ----------------------------------------------------------------------
2463 // ZYAM (for partially negative distributions)
2464 // --> we subtract always the minimum value
2465 Double_t zeroYield = 0.;
2466 Double_t zeroYieldCur = -FLT_MAX;
2468 for(Int_t iMin = 0; iMin<2; iMin++){
2469 zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2470 zeroYield += zeroYieldCur;
2473 //zeroYield = gHist->GetMinimum();
2475 // ----------------------------------------------------------------------
2476 // first calculate the mean
2478 Double_t fWeightedAverage = 0.;
2479 Double_t fNormalization = 0.;
2481 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2483 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2484 if(fVariable == 2 &&
2485 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2486 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2490 fWeightedAverage += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2491 fNormalization += (gHist->GetBinContent(i)-zeroYield);
2494 mean = fWeightedAverage / fNormalization;
2496 // ----------------------------------------------------------------------
2497 // then calculate the higher moments
2508 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2510 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2511 if(fVariable == 2 &&
2512 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2513 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2517 fMu += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2518 fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2519 fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2520 fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2521 fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2522 fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2523 fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2524 fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
2527 // normalize to bin entries!
2528 fMu /= fNormalization;
2529 fMu2 /= fNormalization;
2530 fMu3 /= fNormalization;
2531 fMu4 /= fNormalization;
2532 fMu5 /= fNormalization;
2533 fMu6 /= fNormalization;
2534 fMu7 /= fNormalization;
2535 fMu8 /= fNormalization;
2537 sigma = TMath::Sqrt(fMu2);
2538 skewness = fMu3 / TMath::Power(sigma,3);
2539 kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
2541 // normalize with sigma^r
2542 fMu /= TMath::Power(sigma,1);
2543 fMu2 /= TMath::Power(sigma,2);
2544 fMu3 /= TMath::Power(sigma,3);
2545 fMu4 /= TMath::Power(sigma,4);
2546 fMu5 /= TMath::Power(sigma,5);
2547 fMu6 /= TMath::Power(sigma,6);
2548 fMu7 /= TMath::Power(sigma,7);
2549 fMu8 /= TMath::Power(sigma,8);
2551 // ----------------------------------------------------------------------
2552 // then calculate the higher moment errors
2553 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2554 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2555 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
2557 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2559 // // assuming normal distribution (as done in ROOT)
2560 // meanError = sigma / TMath::Sqrt(normError);
2561 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2562 // skewnessError = TMath::Sqrt(6./(normError));
2563 // kurtosisError = TMath::Sqrt(24./(normError));
2565 // use delta theorem paper (Luo - arXiv:1109.0593v1)
2566 Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2567 Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2568 Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
2569 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2570 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2571 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
2573 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2574 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
2576 if (TMath::Sqrt(normError) != 0){
2577 meanError = sigma / TMath::Sqrt(normError);
2578 sigmaError = TMath::Sqrt(Lambda11);
2579 skewnessError = TMath::Sqrt(Lambda22);
2580 kurtosisError = TMath::Sqrt(Lambda33);
2585 else success = kFALSE;
2592 //____________________________________________________________________//
2593 Float_t AliBalancePsi::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign) {
2595 // calculates dphistar
2597 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2599 static const Double_t kPi = TMath::Pi();
2602 // if (dphistar > 2 * kPi)
2603 // dphistar -= 2 * kPi;
2604 // if (dphistar < -2 * kPi)
2605 // dphistar += 2 * kPi;
2608 dphistar = kPi * 2 - dphistar;
2609 if (dphistar < -kPi)
2610 dphistar = -kPi * 2 - dphistar;
2611 if (dphistar > kPi) // might look funny but is needed
2612 dphistar = kPi * 2 - dphistar;
2617 //____________________________________________________________________//
2618 Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2620 // This method is a copy from AliUEHist::GetBinning
2621 // takes the binning from <configuration> identified by <tag>
2622 // configuration syntax example:
2623 // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
2626 // returns bin edges which have to be deleted by the caller
2628 TString config(configuration);
2629 TObjArray* lines = config.Tokenize("\n");
2630 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2632 TString line(lines->At(i)->GetName());
2633 if (line.BeginsWith(TString(tag) + ":"))
2635 line.Remove(0, strlen(tag) + 1);
2636 line.ReplaceAll(" ", "");
2637 TObjArray* binning = line.Tokenize(",");
2638 Double_t* bins = new Double_t[binning->GetEntriesFast()];
2639 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2640 bins[j] = TString(binning->At(j)->GetName()).Atof();
2642 nBins = binning->GetEntriesFast() - 1;
2651 AliFatal(Form("Tag %s not found in %s", tag, configuration));