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>
37 #include "AliVParticle.h"
38 #include "AliMCParticle.h"
39 #include "AliESDtrack.h"
40 #include "AliAODTrack.h"
42 #include "AliAnalysisTaskTriggeredBF.h"
44 #include "AliBalancePsi.h"
49 ClassImp(AliBalancePsi)
51 //____________________________________________________________________//
52 AliBalancePsi::AliBalancePsi() :
55 fAnalysisLevel("ESD"),
68 fHistConversionbefore(0),
69 fHistConversionafter(0),
71 fHistResonancesBefore(0),
72 fHistResonancesRho(0),
74 fHistResonancesLambda(0),
79 fResonancesCut(kFALSE),
82 fConversionCut(kFALSE),
83 fInvMassCutConversion(0.04),
86 fVertexBinning(kFALSE),
89 fEventClass("EventPlane"){
90 // Default constructor
93 //____________________________________________________________________//
94 AliBalancePsi::AliBalancePsi(const AliBalancePsi& balance):
95 TObject(balance), fShuffle(balance.fShuffle),
96 fAnalysisLevel(balance.fAnalysisLevel),
97 fAnalyzedEvents(balance.fAnalyzedEvents),
98 fCentralityId(balance.fCentralityId),
99 fCentStart(balance.fCentStart),
100 fCentStop(balance.fCentStop),
101 fHistP(balance.fHistP),
102 fHistN(balance.fHistN),
103 fHistPN(balance.fHistPN),
104 fHistNP(balance.fHistNP),
105 fHistPP(balance.fHistPP),
106 fHistNN(balance.fHistNN),
107 fHistHBTbefore(balance.fHistHBTbefore),
108 fHistHBTafter(balance.fHistHBTafter),
109 fHistConversionbefore(balance.fHistConversionbefore),
110 fHistConversionafter(balance.fHistConversionafter),
111 fHistPsiMinusPhi(balance.fHistPsiMinusPhi),
112 fHistResonancesBefore(balance.fHistResonancesBefore),
113 fHistResonancesRho(balance.fHistResonancesRho),
114 fHistResonancesK0(balance.fHistResonancesK0),
115 fHistResonancesLambda(balance.fHistResonancesLambda),
116 fHistQbefore(balance.fHistQbefore),
117 fHistQafter(balance.fHistQafter),
118 fPsiInterval(balance.fPsiInterval),
119 fDeltaEtaMax(balance.fDeltaEtaMax),
120 fResonancesCut(balance.fResonancesCut),
121 fHBTCut(balance.fHBTCut),
122 fHBTCutValue(balance.fHBTCutValue),
123 fConversionCut(balance.fConversionCut),
124 fInvMassCutConversion(balance.fInvMassCutConversion),
125 fQCut(balance.fQCut),
126 fDeltaPtMin(balance.fDeltaPtMin),
127 fVertexBinning(balance.fVertexBinning),
128 fCustomBinning(balance.fCustomBinning),
129 fBinningString(balance.fBinningString),
130 fEventClass("EventPlane"){
134 //____________________________________________________________________//
135 AliBalancePsi::~AliBalancePsi() {
144 delete fHistHBTbefore;
145 delete fHistHBTafter;
146 delete fHistConversionbefore;
147 delete fHistConversionafter;
148 delete fHistPsiMinusPhi;
149 delete fHistResonancesBefore;
150 delete fHistResonancesRho;
151 delete fHistResonancesK0;
152 delete fHistResonancesLambda;
158 //____________________________________________________________________//
159 void AliBalancePsi::InitHistograms() {
160 // single particle histograms
162 // global switch disabling the reference
163 // (to avoid "Replacing existing TH1" if several wagons are created in train)
164 Bool_t oldStatus = TH1::AddDirectoryStatus();
165 TH1::AddDirectory(kFALSE);
167 Int_t anaSteps = 1; // analysis steps
168 Int_t iBinSingle[kTrackVariablesSingle]; // binning for track variables
169 Double_t* dBinsSingle[kTrackVariablesSingle]; // bins for track variables
170 TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
172 // two particle histograms
173 Int_t iBinPair[kTrackVariablesPair]; // binning for track variables
174 Double_t* dBinsPair[kTrackVariablesPair]; // bins for track variables
175 TString axisTitlePair[kTrackVariablesPair]; // axis titles for track variables
179 // =========================================================
180 // The default string (from older versions of AliBalancePsi)
181 // =========================================================
182 TString defaultBinningStr;
183 defaultBinningStr = "multiplicity: 0,10,20,30,40,50,60,70,80,100,100000\n" // Multiplicity Bins
184 "centrality: 0.,5.,10.,20.,30.,40.,50.,60.,70.,80.\n" // Centrality Bins
185 "centralityVertex: 0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.\n" // Centrality Bins (Vertex Binning)
186 "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))
187 "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
188 "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)
189 "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
190 "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
191 "pTVertex: 0.2,1.0,2.0,3.0,4.0,8.0,15.0\n" // pT Bins (Vertex Binning)
192 "vertex: -10., 10.\n" // Vertex Bins
193 "vertexVertex: -10., -7., -5., -3., -1., 1., 3., 5., 7., 10.\n" // Vertex Bins (Vertex Binning)
197 // =========================================================
198 // Customization (adopted from AliUEHistograms)
199 // =========================================================
201 TObjArray* lines = defaultBinningStr.Tokenize("\n");
202 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
204 TString line(lines->At(i)->GetName());
205 TString tag = line(0, line.Index(":")+1);
206 if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
207 fBinningString += line + "\n";
209 AliInfo(Form("Using custom binning for %s", tag.Data()));
212 fBinningString += fCustomBinning;
214 AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
217 // =========================================================
219 // =========================================================
221 //Depending on fEventClass Variable, do one thing or the other...
222 /**********************************************************
224 ======> Modification: Change Event Classification Scheme
226 ---> fEventClass == "EventPlane"
228 Default operation with Event Plane
230 ---> fEventClass == "Multiplicity"
232 Work with reference multiplicity (from GetReferenceMultiplicity, which one is decided in the configuration)
234 ---> fEventClass == "Centrality"
236 Work with Centrality Bins
238 ***********************************************************/
239 if(fEventClass == "Multiplicity"){
240 dBinsSingle[0] = GetBinning(fBinningString, "multiplicity", iBinSingle[0]);
241 dBinsPair[0] = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
242 axisTitleSingle[0] = "reference multiplicity";
243 axisTitlePair[0] = "reference multiplicity";
245 if(fEventClass == "Centrality"){
246 // fine binning in case of vertex Z binning
248 dBinsSingle[0] = GetBinning(fBinningString, "centralityVertex", iBinSingle[0]);
249 dBinsPair[0] = GetBinning(fBinningString, "centralityVertex", iBinPair[0]);
252 dBinsSingle[0] = GetBinning(fBinningString, "centrality", iBinSingle[0]);
253 dBinsPair[0] = GetBinning(fBinningString, "centrality", iBinPair[0]);
255 axisTitleSingle[0] = "Centrality percentile [%]";
256 axisTitlePair[0] = "Centrality percentile [%]";
258 if(fEventClass == "EventPlane"){
259 dBinsSingle[0] = GetBinning(fBinningString, "eventPlane", iBinSingle[0]);
260 dBinsPair[0] = GetBinning(fBinningString, "eventPlane", iBinPair[0]);
261 axisTitleSingle[0] = "#varphi - #Psi_{2} (a.u.)";
262 axisTitlePair[0] = "#varphi - #Psi_{2} (a.u.)";
266 // Delta Eta and Delta Phi
267 // (coarse binning in case of vertex Z binning)
269 dBinsPair[1] = GetBinning(fBinningString, "deltaEtaVertex", iBinPair[1]);
272 dBinsPair[1] = GetBinning(fBinningString, "deltaEta", iBinPair[1]);
274 axisTitlePair[1] = "#Delta#eta";
276 dBinsPair[2] = GetBinning(fBinningString, "deltaPhi", iBinPair[2]);
277 axisTitlePair[2] = "#Delta#varphi (rad)";
280 // pT Trig and pT Assoc
281 // (coarse binning in case of vertex Z binning)
283 dBinsSingle[1] = GetBinning(fBinningString, "pTVertex", iBinSingle[1]);
284 dBinsPair[3] = GetBinning(fBinningString, "pTVertex", iBinPair[3]);
285 dBinsPair[4] = GetBinning(fBinningString, "pTVertex", iBinPair[4]);
288 dBinsSingle[1] = GetBinning(fBinningString, "pT", iBinSingle[1]);
289 dBinsPair[3] = GetBinning(fBinningString, "pT", iBinPair[3]);
290 dBinsPair[4] = GetBinning(fBinningString, "pT", iBinPair[4]);
293 axisTitleSingle[1] = "p_{T,trig.} (GeV/c)";
294 axisTitlePair[3] = "p_{T,trig.} (GeV/c)";
295 axisTitlePair[4] = "p_{T,assoc.} (GeV/c)";
298 // vertex Z binning or not
300 dBinsSingle[2] = GetBinning(fBinningString, "vertexVertex", iBinSingle[2]);
301 dBinsPair[5] = GetBinning(fBinningString, "vertexVertex", iBinPair[5]);
304 dBinsSingle[2] = GetBinning(fBinningString, "vertex", iBinSingle[2]);
305 dBinsPair[5] = GetBinning(fBinningString, "vertex", iBinPair[5]);
308 axisTitleSingle[2] = "v_{Z} (cm)";
309 axisTitlePair[5] = "v_{Z} (cm)";
313 // =========================================================
314 // Create the Output objects (AliTHn)
315 // =========================================================
318 //+ triggered particles
320 if(fShuffle) histName.Append("_shuffle");
321 if(fCentralityId) histName += fCentralityId.Data();
322 fHistP = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
323 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
324 fHistP->SetBinLimits(j, dBinsSingle[j]);
325 fHistP->SetVarTitle(j, axisTitleSingle[j]);
328 //- triggered particles
330 if(fShuffle) histName.Append("_shuffle");
331 if(fCentralityId) histName += fCentralityId.Data();
332 fHistN = new AliTHn(histName.Data(),histName.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
333 for (Int_t j=0; j<kTrackVariablesSingle; j++) {
334 fHistN->SetBinLimits(j, dBinsSingle[j]);
335 fHistN->SetVarTitle(j, axisTitleSingle[j]);
339 histName = "fHistPN";
340 if(fShuffle) histName.Append("_shuffle");
341 if(fCentralityId) histName += fCentralityId.Data();
342 fHistPN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
343 for (Int_t j=0; j<kTrackVariablesPair; j++) {
344 fHistPN->SetBinLimits(j, dBinsPair[j]);
345 fHistPN->SetVarTitle(j, axisTitlePair[j]);
349 histName = "fHistNP";
350 if(fShuffle) histName.Append("_shuffle");
351 if(fCentralityId) histName += fCentralityId.Data();
352 fHistNP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
353 for (Int_t j=0; j<kTrackVariablesPair; j++) {
354 fHistNP->SetBinLimits(j, dBinsPair[j]);
355 fHistNP->SetVarTitle(j, axisTitlePair[j]);
359 histName = "fHistPP";
360 if(fShuffle) histName.Append("_shuffle");
361 if(fCentralityId) histName += fCentralityId.Data();
362 fHistPP = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
363 for (Int_t j=0; j<kTrackVariablesPair; j++) {
364 fHistPP->SetBinLimits(j, dBinsPair[j]);
365 fHistPP->SetVarTitle(j, axisTitlePair[j]);
369 histName = "fHistNN";
370 if(fShuffle) histName.Append("_shuffle");
371 if(fCentralityId) histName += fCentralityId.Data();
372 fHistNN = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
373 for (Int_t j=0; j<kTrackVariablesPair; j++) {
374 fHistNN->SetBinLimits(j, dBinsPair[j]);
375 fHistNN->SetVarTitle(j, axisTitlePair[j]);
378 AliInfo("Finished setting up the AliTHn");
381 fHistHBTbefore = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,2.*TMath::Pi());
382 fHistHBTafter = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,2.*TMath::Pi());
383 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);
384 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);
385 fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
386 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);
387 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);
388 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);
389 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);
390 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);
391 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);
393 TH1::AddDirectory(oldStatus);
397 //____________________________________________________________________//
398 void AliBalancePsi::CalculateBalance(Double_t gReactionPlane,
399 TObjArray *particles,
400 TObjArray *particlesMixed,
402 Double_t kMultorCent,
404 // Calculates the balance function
407 // Initialize histograms if not done yet
409 AliWarning("Histograms not yet initialized! --> Will be done now");
410 AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
414 Double_t trackVariablesSingle[kTrackVariablesSingle];
415 Double_t trackVariablesPair[kTrackVariablesPair];
418 AliWarning("particles TObjArray is NULL pointer --> return");
422 // define end of particle loops
423 Int_t iMax = particles->GetEntriesFast();
426 jMax = particlesMixed->GetEntriesFast();
428 // Eta() is extremely time consuming, therefore cache it for the inner loop here:
429 TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
431 TArrayF secondEta(jMax);
432 TArrayF secondPhi(jMax);
433 TArrayF secondPt(jMax);
434 TArrayS secondCharge(jMax);
435 TArrayD secondCorrection(jMax);
437 for (Int_t i=0; i<jMax; i++){
438 secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
439 secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
440 secondPt[i] = ((AliVParticle*) particlesSecond->At(i))->Pt();
441 secondCharge[i] = (Short_t)((AliVParticle*) particlesSecond->At(i))->Charge();
442 secondCorrection[i] = (Double_t)((AliBFBasicParticle*) particlesSecond->At(i))->Correction(); //==========================correction
445 //TLorenzVector implementation for resonances
446 TLorentzVector vectorMother, vectorDaughter[2];
447 TParticle pPion, pProton, pRho0, pK0s, pLambda;
448 pPion.SetPdgCode(211); //pion
449 pRho0.SetPdgCode(113); //rho0
450 pK0s.SetPdgCode(310); //K0s
451 pProton.SetPdgCode(2212); //proton
452 pLambda.SetPdgCode(3122); //Lambda
453 Double_t gWidthForRho0 = 0.01;
454 Double_t gWidthForK0s = 0.01;
455 Double_t gWidthForLambda = 0.006;
456 Double_t nSigmaRejection = 3.0;
459 for (Int_t i = 0; i < iMax; i++) {
460 //AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
461 AliBFBasicParticle* firstParticle = (AliBFBasicParticle*) particles->At(i); //==========================correction
464 Float_t firstEta = firstParticle->Eta();
465 Float_t firstPhi = firstParticle->Phi();
466 Float_t firstPt = firstParticle->Pt();
467 Float_t firstCorrection = firstParticle->Correction();//==========================correction
469 // Event plane (determine psi bin)
470 Double_t gPsiMinusPhi = 0.;
471 Double_t gPsiMinusPhiBin = -10.;
472 gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
474 if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
475 ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
476 gPsiMinusPhiBin = 0.0;
478 else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
479 ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
480 ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
481 ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
482 gPsiMinusPhiBin = 1.0;
484 else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
485 ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
486 gPsiMinusPhiBin = 2.0;
489 gPsiMinusPhiBin = 3.0;
491 fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
493 Short_t charge1 = (Short_t) firstParticle->Charge();
495 trackVariablesSingle[0] = gPsiMinusPhiBin;
496 trackVariablesSingle[1] = firstPt;
497 if(fEventClass=="Multiplicity" || fEventClass == "Centrality" ) trackVariablesSingle[0] = kMultorCent;
498 trackVariablesSingle[2] = vertexZ;
501 //fill single particle histograms
502 if(charge1 > 0) fHistP->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
503 else if(charge1 < 0) fHistN->Fill(trackVariablesSingle,0,firstCorrection); //==========================correction
506 for(Int_t j = 0; j < jMax; j++) {
508 if(!particlesMixed && j == i) continue; // no auto correlations (only for non mixing)
510 // pT,Assoc < pT,Trig
511 if(firstPt < secondPt[j]) continue;
513 Short_t charge2 = secondCharge[j];
515 trackVariablesPair[0] = trackVariablesSingle[0];
516 trackVariablesPair[1] = firstEta - secondEta[j]; // delta eta
517 trackVariablesPair[2] = firstPhi - secondPhi[j]; // delta phi
518 //if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
519 //trackVariablesPair[2] -= 360.;
520 //if (trackVariablesPair[2] < - 180.)
521 //trackVariablesPair[2] += 360.;
522 if (trackVariablesPair[2] > TMath::Pi()) // delta phi between -pi and pi
523 trackVariablesPair[2] -= 2.*TMath::Pi();
524 if (trackVariablesPair[2] < - TMath::Pi())
525 trackVariablesPair[2] += 2.*TMath::Pi();
526 if (trackVariablesPair[2] < - TMath::Pi()/2.)
527 trackVariablesPair[2] += 2.*TMath::Pi();
529 trackVariablesPair[3] = firstPt; // pt trigger
530 trackVariablesPair[4] = secondPt[j]; // pt
531 trackVariablesPair[5] = vertexZ; // z of the primary vertex
533 //Exclude resonances for the calculation of pairs by looking
534 //at the invariant mass and not considering the pairs that
535 //fall within 3sigma from the mass peak of: rho0, K0s, Lambda
537 if (charge1 * charge2 < 0) {
540 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
541 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
542 vectorMother = vectorDaughter[0] + vectorDaughter[1];
543 fHistResonancesBefore->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
544 if(TMath::Abs(vectorMother.M() - pRho0.GetMass()) <= nSigmaRejection*gWidthForRho0)
546 fHistResonancesRho->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
549 if(TMath::Abs(vectorMother.M() - pK0s.GetMass()) <= nSigmaRejection*gWidthForK0s)
551 fHistResonancesK0->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
555 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pPion.GetMass());
556 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pProton.GetMass());
557 vectorMother = vectorDaughter[0] + vectorDaughter[1];
558 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
561 vectorDaughter[0].SetPtEtaPhiM(firstPt,firstEta,firstPhi,pProton.GetMass());
562 vectorDaughter[1].SetPtEtaPhiM(secondPt[j],secondEta[j],secondPhi[j],pPion.GetMass());
563 vectorMother = vectorDaughter[0] + vectorDaughter[1];
564 if(TMath::Abs(vectorMother.M() - pLambda.GetMass()) <= nSigmaRejection*gWidthForLambda)
566 fHistResonancesLambda->Fill(trackVariablesPair[1],trackVariablesPair[2],vectorMother.M());
572 //if(fHBTCut){ // VERSION 3 (all pairs)
573 if(fHBTCut && charge1 * charge2 > 0){ // VERSION 2 (only for LS)
574 //if( dphi < 3 || deta < 0.01 ){ // VERSION 1
577 Double_t deta = firstEta - secondEta[j];
578 Double_t dphi = firstPhi - secondPhi[j];
579 // VERSION 2 (Taken from DPhiCorrelations)
580 // the variables & cuthave been developed by the HBT group
581 // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
582 fHistHBTbefore->Fill(deta,dphi);
585 if (TMath::Abs(deta) < fHBTCutValue * 2.5 * 3) //fHBTCutValue = 0.02 [default for dphicorrelations]
588 //Float_t phi1rad = firstPhi*TMath::DegToRad();
589 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
590 Float_t phi1rad = firstPhi;
591 Float_t phi2rad = secondPhi[j];
593 // check first boundaries to see if is worth to loop and find the minimum
594 Float_t dphistar1 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 0.8, bSign);
595 Float_t dphistar2 = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, 2.5, bSign);
597 const Float_t kLimit = fHBTCutValue * 3;
599 Float_t dphistarminabs = 1e5;
600 Float_t dphistarmin = 1e5;
602 if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 ) {
603 for (Double_t rad=0.8; rad<2.51; rad+=0.01) {
604 Float_t dphistar = GetDPhiStar(phi1rad, firstPt, charge1, phi2rad, secondPt[j], charge2, rad, bSign);
605 Float_t dphistarabs = TMath::Abs(dphistar);
607 if (dphistarabs < dphistarminabs) {
608 dphistarmin = dphistar;
609 dphistarminabs = dphistarabs;
613 if (dphistarminabs < fHBTCutValue && TMath::Abs(deta) < fHBTCutValue) {
614 //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));
619 fHistHBTafter->Fill(deta,dphi);
624 if (charge1 * charge2 < 0) {
625 Double_t deta = firstEta - secondEta[j];
626 Double_t dphi = firstPhi - secondPhi[j];
628 Float_t m0 = 0.510e-3;
629 Float_t tantheta1 = 1e10;
632 //Float_t phi1rad = firstPhi*TMath::DegToRad();
633 //Float_t phi2rad = secondPhi[j]*TMath::DegToRad();
634 Float_t phi1rad = firstPhi;
635 Float_t phi2rad = secondPhi[j];
637 if (firstEta < -1e-10 || firstEta > 1e-10)
638 tantheta1 = 2 * TMath::Exp(-firstEta) / ( 1 - TMath::Exp(-2*firstEta));
640 Float_t tantheta2 = 1e10;
641 if (secondEta[j] < -1e-10 || secondEta[j] > 1e-10)
642 tantheta2 = 2 * TMath::Exp(-secondEta[j]) / ( 1 - TMath::Exp(-2*secondEta[j]));
644 Float_t e1squ = m0 * m0 + firstPt * firstPt * (1.0 + 1.0 / tantheta1 / tantheta1);
645 Float_t e2squ = m0 * m0 + secondPt[j] * secondPt[j] * (1.0 + 1.0 / tantheta2 / tantheta2);
647 Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( firstPt * secondPt[j] * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );
649 fHistConversionbefore->Fill(deta,dphi,masssqu);
651 if (masssqu < fInvMassCutConversion*fInvMassCutConversion){
652 //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));
655 fHistConversionafter->Fill(deta,dphi,masssqu);
659 // momentum difference cut - suppress femtoscopic effects
662 //Double_t ptMin = 0.1; //const for the time being (should be changeable later on)
663 Double_t ptDifference = TMath::Abs( firstPt - secondPt[j]);
665 fHistQbefore->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
666 if(ptDifference < fDeltaPtMin) continue;
667 fHistQafter->Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
671 if( charge1 > 0 && charge2 < 0) fHistPN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]); //==========================correction
672 else if( charge1 < 0 && charge2 > 0) fHistNP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
673 else if( charge1 > 0 && charge2 > 0) fHistPP->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
674 else if( charge1 < 0 && charge2 < 0) fHistNN->Fill(trackVariablesPair,0,firstCorrection*secondCorrection[j]);//==========================correction
676 //AliWarning(Form("Wrong charge combination: charge1 = %d and charge2 = %d",charge,charge2));
679 }//end of 2nd particle loop
680 }//end of 1st particle loop
683 //____________________________________________________________________//
684 TH1D *AliBalancePsi::GetBalanceFunctionHistogram(Int_t iVariableSingle,
690 Double_t ptTriggerMin,
691 Double_t ptTriggerMax,
692 Double_t ptAssociatedMin,
693 Double_t ptAssociatedMax) {
694 //Returns the BF histogram, extracted from the 6 AliTHn objects
695 //(private members) of the AliBalancePsi class.
696 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
697 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
700 if(psiMin > psiMax-0.00001){
701 AliError("psiMax <= psiMin");
704 if(vertexZMin > vertexZMax-0.00001){
705 AliError("vertexZMax <= vertexZMin");
708 if(ptTriggerMin > ptTriggerMax-0.00001){
709 AliError("ptTriggerMax <= ptTriggerMin");
712 if(ptAssociatedMin > ptAssociatedMax-0.00001){
713 AliError("ptAssociatedMax <= ptAssociatedMin");
718 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
719 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
720 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
721 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
722 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
723 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
727 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
728 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
729 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
730 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
731 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
732 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
736 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
737 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
738 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
739 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
740 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
741 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
742 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
746 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
747 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
748 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
749 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
750 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
753 //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));
755 // Project into the wanted space (1st: analysis step, 2nd: axis)
756 TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,iVariablePair); //
757 TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,iVariablePair); //
758 TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,iVariablePair); //
759 TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,iVariablePair); //
760 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle); //
761 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle); //
763 TH1D *gHistBalanceFunctionHistogram = 0x0;
764 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
765 gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
766 gHistBalanceFunctionHistogram->Reset();
768 switch(iVariablePair) {
770 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
771 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
774 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
775 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
785 hTemp1->Add(hTemp3,-1.);
786 hTemp1->Scale(1./hTemp5->Integral());
787 hTemp2->Add(hTemp4,-1.);
788 hTemp2->Scale(1./hTemp6->Integral());
789 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
790 gHistBalanceFunctionHistogram->Scale(0.5);
792 //normalize to bin width
793 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
796 return gHistBalanceFunctionHistogram;
799 //____________________________________________________________________//
800 TH1D *AliBalancePsi::GetBalanceFunctionHistogram2pMethod(Int_t iVariableSingle,
806 Double_t ptTriggerMin,
807 Double_t ptTriggerMax,
808 Double_t ptAssociatedMin,
809 Double_t ptAssociatedMax,
810 AliBalancePsi *bfMix) {
811 //Returns the BF histogram, extracted from the 6 AliTHn objects
812 //after dividing each correlation function by the Event Mixing one
813 //(private members) of the AliBalancePsi class.
814 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
815 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
818 if(psiMin > psiMax-0.00001){
819 AliError("psiMax <= psiMin");
822 if(vertexZMin > vertexZMax-0.00001){
823 AliError("vertexZMax <= vertexZMin");
826 if(ptTriggerMin > ptTriggerMax-0.00001){
827 AliError("ptTriggerMax <= ptTriggerMin");
830 if(ptAssociatedMin > ptAssociatedMax-0.00001){
831 AliError("ptAssociatedMax <= ptAssociatedMin");
835 // pt trigger (fixed for all vertices/psi/centralities)
836 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
837 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
838 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
839 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
840 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
841 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
842 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
845 // pt associated (fixed for all vertices/psi/centralities)
846 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
847 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
848 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
849 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
850 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
853 // ============================================================================================
854 // the same for event mixing
855 AliTHn *fHistPMix = bfMix->GetHistNp();
856 AliTHn *fHistNMix = bfMix->GetHistNn();
857 AliTHn *fHistPNMix = bfMix->GetHistNpn();
858 AliTHn *fHistNPMix = bfMix->GetHistNnp();
859 AliTHn *fHistPPMix = bfMix->GetHistNpp();
860 AliTHn *fHistNNMix = bfMix->GetHistNnn();
862 // pt trigger (fixed for all vertices/psi/centralities)
863 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
864 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
865 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
866 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
867 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
868 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
869 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
872 // pt associated (fixed for all vertices/psi/centralities)
873 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
874 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
875 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
876 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
877 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
880 // ============================================================================================
881 // ranges (in bins) for vertexZ and centrality (psi)
882 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
883 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
884 Int_t binVertexMin = 0;
885 Int_t binVertexMax = 0;
887 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
888 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
890 Double_t binPsiLowEdge = 0.;
891 Double_t binPsiUpEdge = 1000.;
892 Double_t binVertexLowEdge = 0.;
893 Double_t binVertexUpEdge = 1000.;
900 // loop over all bins
901 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
902 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
904 cout<<"In the balance function (1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
906 // determine the bin edges for this bin
907 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
908 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
910 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
911 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
915 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
916 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
917 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
918 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
919 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
920 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
924 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
925 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
926 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
927 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
928 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
929 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
932 // ============================================================================================
933 // the same for event mixing
936 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
937 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
938 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
939 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
940 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
941 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
945 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
946 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
947 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
948 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
949 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
950 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
953 // ============================================================================================
955 //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));
957 // Project into the wanted space (1st: analysis step, 2nd: axis)
958 TH1D* hTempHelper1 = (TH1D*)fHistPN->Project(0,iVariablePair);
959 TH1D* hTempHelper2 = (TH1D*)fHistNP->Project(0,iVariablePair);
960 TH1D* hTempHelper3 = (TH1D*)fHistPP->Project(0,iVariablePair);
961 TH1D* hTempHelper4 = (TH1D*)fHistNN->Project(0,iVariablePair);
962 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,iVariableSingle);
963 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,iVariableSingle);
965 // ============================================================================================
966 // the same for event mixing
967 TH1D* hTempHelper1Mix = (TH1D*)fHistPNMix->Project(0,iVariablePair);
968 TH1D* hTempHelper2Mix = (TH1D*)fHistNPMix->Project(0,iVariablePair);
969 TH1D* hTempHelper3Mix = (TH1D*)fHistPPMix->Project(0,iVariablePair);
970 TH1D* hTempHelper4Mix = (TH1D*)fHistNNMix->Project(0,iVariablePair);
971 TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,iVariableSingle);
972 TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,iVariableSingle);
973 // ============================================================================================
975 hTempHelper1->Sumw2();
976 hTempHelper2->Sumw2();
977 hTempHelper3->Sumw2();
978 hTempHelper4->Sumw2();
981 hTempHelper1Mix->Sumw2();
982 hTempHelper2Mix->Sumw2();
983 hTempHelper3Mix->Sumw2();
984 hTempHelper4Mix->Sumw2();
988 // first put everything on positive x - axis (to be comparable with published data) --> ONLY IN THIS OPTION!
990 Double_t helperEndBin = 1.6;
991 if(iVariablePair==2) helperEndBin = TMath::Pi();
993 TH1D* hTempPos1 = new TH1D(Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
994 TH1D* hTempPos2 = new TH1D(Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
995 TH1D* hTempPos3 = new TH1D(Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
996 TH1D* hTempPos4 = new TH1D(Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
997 TH1D* hTempPos1Mix = new TH1D(Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
998 TH1D* hTempPos2Mix = new TH1D(Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
999 TH1D* hTempPos3Mix = new TH1D(Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
1000 TH1D* hTempPos4Mix = new TH1D(Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTempPos4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1002 TH1D* hTemp1 = new TH1D(Form("hTemp1_%d_%d",iBinPsi,iBinVertex),Form("hTemp1_%d_%d",iBinPsi,iBinVertex),hTempHelper1->GetNbinsX()/2,0,helperEndBin);
1003 TH1D* hTemp2 = new TH1D(Form("hTemp2_%d_%d",iBinPsi,iBinVertex),Form("hTemp2_%d_%d",iBinPsi,iBinVertex),hTempHelper2->GetNbinsX()/2,0,helperEndBin);
1004 TH1D* hTemp3 = new TH1D(Form("hTemp3_%d_%d",iBinPsi,iBinVertex),Form("hTemp3_%d_%d",iBinPsi,iBinVertex),hTempHelper3->GetNbinsX()/2,0,helperEndBin);
1005 TH1D* hTemp4 = new TH1D(Form("hTemp4_%d_%d",iBinPsi,iBinVertex),Form("hTemp4_%d_%d",iBinPsi,iBinVertex),hTempHelper4->GetNbinsX()/2,0,helperEndBin);
1006 TH1D* hTemp1Mix = new TH1D(Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp1Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper1Mix->GetNbinsX()/2,0,helperEndBin);
1007 TH1D* hTemp2Mix = new TH1D(Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp2Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper2Mix->GetNbinsX()/2,0,helperEndBin);
1008 TH1D* hTemp3Mix = new TH1D(Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp3Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper3Mix->GetNbinsX()/2,0,helperEndBin);
1009 TH1D* hTemp4Mix = new TH1D(Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),Form("hTemp4Mix_%d_%d",iBinPsi,iBinVertex),hTempHelper4Mix->GetNbinsX()/2,0,helperEndBin);
1016 hTempPos1Mix->Sumw2();
1017 hTempPos2Mix->Sumw2();
1018 hTempPos3Mix->Sumw2();
1019 hTempPos4Mix->Sumw2();
1031 for(Int_t i=0;i<hTempHelper1->GetNbinsX();i++){
1033 Double_t binCenter = hTempHelper1->GetXaxis()->GetBinCenter(i+1);
1035 if(iVariablePair==1){
1037 hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1038 hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1039 hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1040 hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1041 hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1042 hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1043 hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1044 hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1046 hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1047 hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1048 hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1049 hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1050 hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1051 hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1052 hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1053 hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1056 hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1057 hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1058 hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1059 hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1060 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1061 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1062 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1063 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1065 hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1066 hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1067 hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1068 hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1069 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1070 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1071 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1072 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1075 else if(iVariablePair==2){
1076 if(binCenter>0 && binCenter<TMath::Pi()){
1077 hTempPos1->SetBinContent(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinContent(i+1));
1078 hTempPos2->SetBinContent(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinContent(i+1));
1079 hTempPos3->SetBinContent(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinContent(i+1));
1080 hTempPos4->SetBinContent(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinContent(i+1));
1081 hTempPos1Mix->SetBinContent(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinContent(i+1));
1082 hTempPos2Mix->SetBinContent(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinContent(i+1));
1083 hTempPos3Mix->SetBinContent(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinContent(i+1));
1084 hTempPos4Mix->SetBinContent(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinContent(i+1));
1086 hTempPos1->SetBinError(hTempPos1->FindBin(binCenter),hTempHelper1->GetBinError(i+1));
1087 hTempPos2->SetBinError(hTempPos2->FindBin(binCenter),hTempHelper2->GetBinError(i+1));
1088 hTempPos3->SetBinError(hTempPos3->FindBin(binCenter),hTempHelper3->GetBinError(i+1));
1089 hTempPos4->SetBinError(hTempPos4->FindBin(binCenter),hTempHelper4->GetBinError(i+1));
1090 hTempPos1Mix->SetBinError(hTempPos1Mix->FindBin(binCenter),hTempHelper1Mix->GetBinError(i+1));
1091 hTempPos2Mix->SetBinError(hTempPos2Mix->FindBin(binCenter),hTempHelper2Mix->GetBinError(i+1));
1092 hTempPos3Mix->SetBinError(hTempPos3Mix->FindBin(binCenter),hTempHelper3Mix->GetBinError(i+1));
1093 hTempPos4Mix->SetBinError(hTempPos4Mix->FindBin(binCenter),hTempHelper4Mix->GetBinError(i+1));
1095 else if(binCenter<0){
1096 hTemp1->SetBinContent(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinContent(i+1));
1097 hTemp2->SetBinContent(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinContent(i+1));
1098 hTemp3->SetBinContent(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinContent(i+1));
1099 hTemp4->SetBinContent(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinContent(i+1));
1100 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1101 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1102 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1103 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1105 hTemp1->SetBinError(hTemp1->FindBin(-binCenter),hTempHelper1->GetBinError(i+1));
1106 hTemp2->SetBinError(hTemp2->FindBin(-binCenter),hTempHelper2->GetBinError(i+1));
1107 hTemp3->SetBinError(hTemp3->FindBin(-binCenter),hTempHelper3->GetBinError(i+1));
1108 hTemp4->SetBinError(hTemp4->FindBin(-binCenter),hTempHelper4->GetBinError(i+1));
1109 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(-binCenter),hTempHelper1Mix->GetBinError(i+1));
1110 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(-binCenter),hTempHelper2Mix->GetBinError(i+1));
1111 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(-binCenter),hTempHelper3Mix->GetBinError(i+1));
1112 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(-binCenter),hTempHelper4Mix->GetBinError(i+1));
1115 hTemp1->SetBinContent(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinContent(i+1));
1116 hTemp2->SetBinContent(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinContent(i+1));
1117 hTemp3->SetBinContent(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinContent(i+1));
1118 hTemp4->SetBinContent(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinContent(i+1));
1119 hTemp1Mix->SetBinContent(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinContent(i+1));
1120 hTemp2Mix->SetBinContent(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinContent(i+1));
1121 hTemp3Mix->SetBinContent(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinContent(i+1));
1122 hTemp4Mix->SetBinContent(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinContent(i+1));
1124 hTemp1->SetBinError(hTemp1->FindBin(TMath::Pi()-binCenter),hTempHelper1->GetBinError(i+1));
1125 hTemp2->SetBinError(hTemp2->FindBin(TMath::Pi()-binCenter),hTempHelper2->GetBinError(i+1));
1126 hTemp3->SetBinError(hTemp3->FindBin(TMath::Pi()-binCenter),hTempHelper3->GetBinError(i+1));
1127 hTemp4->SetBinError(hTemp4->FindBin(TMath::Pi()-binCenter),hTempHelper4->GetBinError(i+1));
1128 hTemp1Mix->SetBinError(hTemp1Mix->FindBin(TMath::Pi()-binCenter),hTempHelper1Mix->GetBinError(i+1));
1129 hTemp2Mix->SetBinError(hTemp2Mix->FindBin(TMath::Pi()-binCenter),hTempHelper2Mix->GetBinError(i+1));
1130 hTemp3Mix->SetBinError(hTemp3Mix->FindBin(TMath::Pi()-binCenter),hTempHelper3Mix->GetBinError(i+1));
1131 hTemp4Mix->SetBinError(hTemp4Mix->FindBin(TMath::Pi()-binCenter),hTempHelper4Mix->GetBinError(i+1));
1136 hTemp1->Add(hTempPos1);
1137 hTemp2->Add(hTempPos2);
1138 hTemp3->Add(hTempPos3);
1139 hTemp4->Add(hTempPos4);
1140 hTemp1Mix->Add(hTempPos1Mix);
1141 hTemp2Mix->Add(hTempPos2Mix);
1142 hTemp3Mix->Add(hTempPos3Mix);
1143 hTemp4Mix->Add(hTempPos4Mix);
1145 hTemp1->Scale(1./hTemp5->Integral());
1146 hTemp3->Scale(1./hTemp5->Integral());
1147 hTemp2->Scale(1./hTemp6->Integral());
1148 hTemp4->Scale(1./hTemp6->Integral());
1150 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1151 // does not work here, so normalize also to trigger particles
1152 // --> careful: gives different integrals then as with full 2D method
1153 hTemp1Mix->Scale(1./hTemp5Mix->Integral());
1154 hTemp3Mix->Scale(1./hTemp5Mix->Integral());
1155 hTemp2Mix->Scale(1./hTemp6Mix->Integral());
1156 hTemp4Mix->Scale(1./hTemp6Mix->Integral());
1158 hTemp1->Divide(hTemp1Mix);
1159 hTemp2->Divide(hTemp2Mix);
1160 hTemp3->Divide(hTemp3Mix);
1161 hTemp4->Divide(hTemp4Mix);
1163 // for the first: clone
1164 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1165 h1 = (TH1D*)hTemp1->Clone();
1166 h2 = (TH1D*)hTemp2->Clone();
1167 h3 = (TH1D*)hTemp3->Clone();
1168 h4 = (TH1D*)hTemp4->Clone();
1170 else{ // otherwise: add for averaging
1190 delete hTempPos1Mix;
1191 delete hTempPos2Mix;
1192 delete hTempPos3Mix;
1193 delete hTempPos4Mix;
1197 TH1D *gHistBalanceFunctionHistogram = 0x0;
1198 if((h1)&&(h2)&&(h3)&&(h4)) {
1200 // average over number of bins nbinsVertex * nbinsPsi
1201 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1202 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1203 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1204 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1206 gHistBalanceFunctionHistogram = (TH1D*)h1->Clone();
1207 gHistBalanceFunctionHistogram->Reset();
1209 switch(iVariablePair) {
1211 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1212 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1215 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1216 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1225 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1226 gHistBalanceFunctionHistogram->Scale(0.5);
1228 //normalize to bin width
1229 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1232 return gHistBalanceFunctionHistogram;
1235 //____________________________________________________________________//
1236 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi(Double_t psiMin,
1238 Double_t vertexZMin,
1239 Double_t vertexZMax,
1240 Double_t ptTriggerMin,
1241 Double_t ptTriggerMax,
1242 Double_t ptAssociatedMin,
1243 Double_t ptAssociatedMax) {
1244 //Returns the BF histogram in Delta eta vs Delta phi,
1245 //extracted from the 6 AliTHn objects
1246 //(private members) of the AliBalancePsi class.
1247 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1248 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1251 if(psiMin > psiMax-0.00001){
1252 AliError("psiMax <= psiMin");
1255 if(vertexZMin > vertexZMax-0.00001){
1256 AliError("vertexZMax <= vertexZMin");
1259 if(ptTriggerMin > ptTriggerMax-0.00001){
1260 AliError("ptTriggerMax <= ptTriggerMin");
1263 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1264 AliError("ptAssociatedMax <= ptAssociatedMin");
1268 TString histName = "gHistBalanceFunctionHistogram2D";
1271 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1272 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1273 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1274 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1275 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1276 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1280 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1281 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1282 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1283 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1284 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1285 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1289 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1290 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1291 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1292 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1293 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1294 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1295 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1299 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1300 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1301 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1302 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1303 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1306 //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)));
1308 // Project into the wanted space (1st: analysis step, 2nd: axis)
1309 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1310 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1311 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1312 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1313 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1314 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1316 TH2D *gHistBalanceFunctionHistogram = 0x0;
1317 if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
1318 gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
1319 gHistBalanceFunctionHistogram->Reset();
1320 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1321 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1322 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1328 hTemp1->Add(hTemp3,-1.);
1329 hTemp1->Scale(1./hTemp5->Integral());
1330 hTemp2->Add(hTemp4,-1.);
1331 hTemp2->Scale(1./hTemp6->Integral());
1332 gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
1333 gHistBalanceFunctionHistogram->Scale(0.5);
1335 //normalize to bin width
1336 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1339 return gHistBalanceFunctionHistogram;
1342 //____________________________________________________________________//
1343 TH2D *AliBalancePsi::GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(Double_t psiMin,
1345 Double_t vertexZMin,
1346 Double_t vertexZMax,
1347 Double_t ptTriggerMin,
1348 Double_t ptTriggerMax,
1349 Double_t ptAssociatedMin,
1350 Double_t ptAssociatedMax,
1351 AliBalancePsi *bfMix) {
1352 //Returns the BF histogram in Delta eta vs Delta phi,
1353 //after dividing each correlation function by the Event Mixing one
1354 //extracted from the 6 AliTHn objects
1355 //(private members) of the AliBalancePsi class.
1356 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1357 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1358 TString histName = "gHistBalanceFunctionHistogram2D";
1361 AliError("balance function object for event mixing not available");
1366 if(psiMin > psiMax-0.00001){
1367 AliError("psiMax <= psiMin");
1370 if(vertexZMin > vertexZMax-0.00001){
1371 AliError("vertexZMax <= vertexZMin");
1374 if(ptTriggerMin > ptTriggerMax-0.00001){
1375 AliError("ptTriggerMax <= ptTriggerMin");
1378 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1379 AliError("ptAssociatedMax <= ptAssociatedMin");
1383 // pt trigger (fixed for all vertices/psi/centralities)
1384 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1385 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1386 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1387 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1388 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1389 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1390 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1393 // pt associated (fixed for all vertices/psi/centralities)
1394 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1395 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1396 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1397 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1398 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1401 // ============================================================================================
1402 // the same for event mixing
1403 AliTHn *fHistPMix = bfMix->GetHistNp();
1404 AliTHn *fHistNMix = bfMix->GetHistNn();
1405 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1406 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1407 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1408 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1410 // pt trigger (fixed for all vertices/psi/centralities)
1411 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1412 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1413 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1414 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1415 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1416 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1417 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1420 // pt associated (fixed for all vertices/psi/centralities)
1421 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1422 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1423 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1424 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1425 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1428 // ============================================================================================
1429 // ranges (in bins) for vertexZ and centrality (psi)
1430 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1431 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1432 Int_t binVertexMin = 0;
1433 Int_t binVertexMax = 0;
1435 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1436 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1438 Double_t binPsiLowEdge = 0.;
1439 Double_t binPsiUpEdge = 0.;
1440 Double_t binVertexLowEdge = 0.;
1441 Double_t binVertexUpEdge = 0.;
1448 // loop over all bins
1449 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1450 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1452 cout<<"In the balance function (2D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1454 // determine the bin edges for this bin
1455 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi);
1456 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi);
1458 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex);
1459 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex);
1463 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1464 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1465 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1466 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1467 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1468 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1472 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1473 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1474 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1475 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1476 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1477 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1482 // ============================================================================================
1483 // the same for event mixing
1486 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1487 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1488 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1489 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1490 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1491 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(binPsiLowEdge,binPsiUpEdge-0.00001);
1495 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1496 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1497 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1498 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1499 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1500 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(binVertexLowEdge,binVertexUpEdge-0.00001);
1503 // ============================================================================================
1506 //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)));
1508 // Project into the wanted space (1st: analysis step, 2nd: axis)
1509 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1510 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1511 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1512 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1513 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1514 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1516 // ============================================================================================
1517 // the same for event mixing
1518 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1519 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1520 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1521 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1522 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1523 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1524 // ============================================================================================
1535 hTemp1->Scale(1./hTemp5->Integral());
1536 hTemp3->Scale(1./hTemp5->Integral());
1537 hTemp2->Scale(1./hTemp6->Integral());
1538 hTemp4->Scale(1./hTemp6->Integral());
1540 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1541 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1542 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1543 hTemp1Mix->Scale(1./mixedNorm1);
1544 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1545 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1546 hTemp2Mix->Scale(1./mixedNorm2);
1547 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1548 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1549 hTemp3Mix->Scale(1./mixedNorm3);
1550 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1551 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1552 hTemp4Mix->Scale(1./mixedNorm4);
1554 hTemp1->Divide(hTemp1Mix);
1555 hTemp2->Divide(hTemp2Mix);
1556 hTemp3->Divide(hTemp3Mix);
1557 hTemp4->Divide(hTemp4Mix);
1559 // for the first: clone
1560 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1561 h1 = (TH2D*)hTemp1->Clone();
1562 h2 = (TH2D*)hTemp2->Clone();
1563 h3 = (TH2D*)hTemp3->Clone();
1564 h4 = (TH2D*)hTemp4->Clone();
1566 else{ // otherwise: add for averaging
1575 TH2D *gHistBalanceFunctionHistogram = 0x0;
1576 if((h1)&&(h2)&&(h3)&&(h4)) {
1578 // average over number of bins nbinsVertex * nbinsPsi
1579 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1580 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1581 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1582 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1584 gHistBalanceFunctionHistogram = (TH2D*)h1->Clone();
1585 gHistBalanceFunctionHistogram->Reset();
1586 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1587 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1588 gHistBalanceFunctionHistogram->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1593 gHistBalanceFunctionHistogram->Add(h1,h2,1.,1.);
1594 gHistBalanceFunctionHistogram->Scale(0.5);
1596 //normalize to bin width
1597 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)*(Double_t)gHistBalanceFunctionHistogram->GetYaxis()->GetBinWidth(1)));
1600 return gHistBalanceFunctionHistogram;
1603 //____________________________________________________________________//
1604 TH1D *AliBalancePsi::GetBalanceFunction1DFrom2D2pMethod(Bool_t bPhi,
1607 Double_t vertexZMin,
1608 Double_t vertexZMax,
1609 Double_t ptTriggerMin,
1610 Double_t ptTriggerMax,
1611 Double_t ptAssociatedMin,
1612 Double_t ptAssociatedMax,
1613 AliBalancePsi *bfMix) {
1614 //Returns the BF histogram in Delta eta OR Delta phi,
1615 //after dividing each correlation function by the Event Mixing one
1616 // (But the division is done here in 2D, this was basically done to check the Integral)
1617 //extracted from the 6 AliTHn objects
1618 //(private members) of the AliBalancePsi class.
1619 //iVariableSingle: 0(phi-Psi), 1(pt-trigger)
1620 //iVariablePair: 0(phi-Psi) 1(Delta eta), 2(Delta phi), 3(pt-trigger), 4(pt-associated
1621 TString histName = "gHistBalanceFunctionHistogram1D";
1624 AliError("balance function object for event mixing not available");
1629 if(psiMin > psiMax-0.00001){
1630 AliError("psiMax <= psiMin");
1633 if(vertexZMin > vertexZMax-0.00001){
1634 AliError("vertexZMax <= vertexZMin");
1637 if(ptTriggerMin > ptTriggerMax-0.00001){
1638 AliError("ptTriggerMax <= ptTriggerMin");
1641 if(ptAssociatedMin > ptAssociatedMax-0.00001){
1642 AliError("ptAssociatedMax <= ptAssociatedMin");
1646 // pt trigger (fixed for all vertices/psi/centralities)
1647 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1648 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1649 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1650 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1651 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1652 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1653 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1656 // pt associated (fixed for all vertices/psi/centralities)
1657 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1658 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1659 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1660 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1661 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1664 // ============================================================================================
1665 // the same for event mixing
1666 AliTHn *fHistPMix = bfMix->GetHistNp();
1667 AliTHn *fHistNMix = bfMix->GetHistNn();
1668 AliTHn *fHistPNMix = bfMix->GetHistNpn();
1669 AliTHn *fHistNPMix = bfMix->GetHistNnp();
1670 AliTHn *fHistPPMix = bfMix->GetHistNpp();
1671 AliTHn *fHistNNMix = bfMix->GetHistNnn();
1673 // pt trigger (fixed for all vertices/psi/centralities)
1674 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
1675 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1676 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1677 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1678 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1679 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1680 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
1683 // pt associated (fixed for all vertices/psi/centralities)
1684 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)) {
1685 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1686 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1687 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1688 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
1691 // ============================================================================================
1692 // ranges (in bins) for vertexZ and centrality (psi)
1693 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin);
1694 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1695 Int_t binVertexMin = 0;
1696 Int_t binVertexMax = 0;
1698 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin);
1699 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1707 // loop over all bins
1708 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1709 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1711 cout<<"In the balance function (2D->1D) loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1714 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1715 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1716 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1717 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1718 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1719 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1723 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1724 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1725 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1726 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1727 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1728 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1731 // ============================================================================================
1732 // the same for event mixing
1735 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1736 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1737 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1738 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1739 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1740 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
1744 fHistPMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1745 fHistNMix->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1746 fHistPNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1747 fHistNPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1748 fHistPPMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1749 fHistNNMix->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
1751 // ============================================================================================
1753 //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)));
1755 // Project into the wanted space (1st: analysis step, 2nd: axis)
1756 TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,1,2);
1757 TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,1,2);
1758 TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,1,2);
1759 TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,1,2);
1760 TH1D* hTemp5 = (TH1D*)fHistP->Project(0,1);
1761 TH1D* hTemp6 = (TH1D*)fHistN->Project(0,1);
1763 // ============================================================================================
1764 // the same for event mixing
1765 TH2D* hTemp1Mix = (TH2D*)fHistPNMix->Project(0,1,2);
1766 TH2D* hTemp2Mix = (TH2D*)fHistNPMix->Project(0,1,2);
1767 TH2D* hTemp3Mix = (TH2D*)fHistPPMix->Project(0,1,2);
1768 TH2D* hTemp4Mix = (TH2D*)fHistNNMix->Project(0,1,2);
1769 // TH1D* hTemp5Mix = (TH1D*)fHistPMix->Project(0,1);
1770 // TH1D* hTemp6Mix = (TH1D*)fHistNMix->Project(0,1);
1771 // ============================================================================================
1782 hTemp1->Scale(1./hTemp5->Integral());
1783 hTemp3->Scale(1./hTemp5->Integral());
1784 hTemp2->Scale(1./hTemp6->Integral());
1785 hTemp4->Scale(1./hTemp6->Integral());
1787 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1788 Double_t mixedNorm1 = hTemp1Mix->Integral(hTemp1Mix->GetXaxis()->FindBin(0-10e-5),hTemp1Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp1Mix->GetNbinsX());
1789 mixedNorm1 /= hTemp1Mix->GetNbinsY()*(hTemp1Mix->GetXaxis()->FindBin(0.01) - hTemp1Mix->GetXaxis()->FindBin(-0.01) + 1);
1790 hTemp1Mix->Scale(1./mixedNorm1);
1791 Double_t mixedNorm2 = hTemp2Mix->Integral(hTemp2Mix->GetXaxis()->FindBin(0-10e-5),hTemp2Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp2Mix->GetNbinsX());
1792 mixedNorm2 /= hTemp2Mix->GetNbinsY()*(hTemp2Mix->GetXaxis()->FindBin(0.01) - hTemp2Mix->GetXaxis()->FindBin(-0.01) + 1);
1793 hTemp2Mix->Scale(1./mixedNorm2);
1794 Double_t mixedNorm3 = hTemp3Mix->Integral(hTemp3Mix->GetXaxis()->FindBin(0-10e-5),hTemp3Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp3Mix->GetNbinsX());
1795 mixedNorm3 /= hTemp3Mix->GetNbinsY()*(hTemp3Mix->GetXaxis()->FindBin(0.01) - hTemp3Mix->GetXaxis()->FindBin(-0.01) + 1);
1796 hTemp3Mix->Scale(1./mixedNorm3);
1797 Double_t mixedNorm4 = hTemp4Mix->Integral(hTemp4Mix->GetXaxis()->FindBin(0-10e-5),hTemp4Mix->GetXaxis()->FindBin(0+10e-5),1,hTemp4Mix->GetNbinsX());
1798 mixedNorm4 /= hTemp4Mix->GetNbinsY()*(hTemp4Mix->GetXaxis()->FindBin(0.01) - hTemp4Mix->GetXaxis()->FindBin(-0.01) + 1);
1799 hTemp4Mix->Scale(1./mixedNorm4);
1801 hTemp1->Divide(hTemp1Mix);
1802 hTemp2->Divide(hTemp2Mix);
1803 hTemp3->Divide(hTemp3Mix);
1804 hTemp4->Divide(hTemp4Mix);
1806 // for the first: clone
1807 if(iBinPsi == binPsiMin && iBinVertex == binVertexMin ){
1808 h1 = (TH2D*)hTemp1->Clone();
1809 h2 = (TH2D*)hTemp2->Clone();
1810 h3 = (TH2D*)hTemp3->Clone();
1811 h4 = (TH2D*)hTemp4->Clone();
1813 else{ // otherwise: add for averaging
1823 TH1D *gHistBalanceFunctionHistogram = 0x0;
1824 if((h1)&&(h2)&&(h3)&&(h4)) {
1826 // average over number of bins nbinsVertex * nbinsPsi
1827 h1->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1828 h2->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1829 h3->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1830 h4->Scale(1./((Double_t)(binPsiMax-binPsiMin+1)*(binVertexMax-binVertexMin+1)));
1832 // now only project on one axis
1833 TH1D *h1DTemp1 = NULL;
1834 TH1D *h1DTemp2 = NULL;
1835 TH1D *h1DTemp3 = NULL;
1836 TH1D *h1DTemp4 = NULL;
1838 h1DTemp1 = (TH1D*)h1->ProjectionX(Form("%s_projX",h1->GetName()));
1839 h1DTemp2 = (TH1D*)h2->ProjectionX(Form("%s_projX",h2->GetName()));
1840 h1DTemp3 = (TH1D*)h3->ProjectionX(Form("%s_projX",h3->GetName()));
1841 h1DTemp4 = (TH1D*)h4->ProjectionX(Form("%s_projX",h4->GetName()));
1844 h1DTemp1 = (TH1D*)h1->ProjectionY(Form("%s_projY",h1->GetName()));
1845 h1DTemp2 = (TH1D*)h2->ProjectionY(Form("%s_projY",h2->GetName()));
1846 h1DTemp3 = (TH1D*)h3->ProjectionY(Form("%s_projY",h3->GetName()));
1847 h1DTemp4 = (TH1D*)h4->ProjectionY(Form("%s_projY",h4->GetName()));
1850 gHistBalanceFunctionHistogram = (TH1D*)h1DTemp1->Clone(Form("%s_clone",h1DTemp1->GetName()));
1851 gHistBalanceFunctionHistogram->Reset();
1853 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#eta");
1854 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#eta)");
1857 gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1858 gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta#varphi)");
1861 h1DTemp1->Add(h1DTemp3,-1.);
1862 h1DTemp2->Add(h1DTemp4,-1.);
1864 gHistBalanceFunctionHistogram->Add(h1DTemp1,h1DTemp2,1.,1.);
1865 gHistBalanceFunctionHistogram->Scale(0.5);
1867 //normalize to bin width
1868 gHistBalanceFunctionHistogram->Scale(1./((Double_t)gHistBalanceFunctionHistogram->GetXaxis()->GetBinWidth(1)));
1871 return gHistBalanceFunctionHistogram;
1875 //____________________________________________________________________//
1876 TH2D *AliBalancePsi::GetCorrelationFunction(TString type,
1879 Double_t vertexZMin,
1880 Double_t vertexZMax,
1881 Double_t ptTriggerMin,
1882 Double_t ptTriggerMax,
1883 Double_t ptAssociatedMin,
1884 Double_t ptAssociatedMax,
1885 AliBalancePsi *bMixed,
1887 Double_t normalizationRangePhi) {
1889 // Returns the 2D correlation function for "type"(PN,NP,PP,NN) pairs,
1890 // does the division by event mixing inside,
1891 // and averaging over several vertexZ and centrality bins
1894 if(type != "PN" && type != "NP" && type != "PP" && type != "NN" && type != "ALL"){
1895 AliError("Only types allowed: PN,NP,PP,NN,ALL");
1899 AliError("No Event Mixing AliTHn");
1905 TH2D *fMixed = NULL;
1907 // ranges (in bins) for vertexZ and centrality (psi)
1908 Int_t binPsiMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMin+0.00001);
1909 Int_t binPsiMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->FindBin(psiMax-0.00001);
1910 Int_t binVertexMin = 0;
1911 Int_t binVertexMax = 0;
1913 binVertexMin = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMin+0.00001);
1914 binVertexMax = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->FindBin(vertexZMax-0.00001);
1916 Double_t binPsiLowEdge = 0.;
1917 Double_t binPsiUpEdge = 1000.;
1918 Double_t binVertexLowEdge = 0.;
1919 Double_t binVertexUpEdge = 1000.;
1921 // if event mixing is empty, we can not add that sub bin
1922 // need to record the number of triggers for correct normalization
1923 Double_t nTrigSubBinEmpty = 0.;
1925 // loop over all bins
1926 for(Int_t iBinPsi = binPsiMin; iBinPsi <= binPsiMax; iBinPsi++){
1927 for(Int_t iBinVertex = binVertexMin; iBinVertex <= binVertexMax; iBinVertex++){
1929 //cout<<"In the correlation function loop: "<<iBinPsi<<" (psiBin), "<<iBinVertex<<" (vertexBin) "<<endl;
1931 // determine the bin edges for this bin
1932 binPsiLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinLowEdge(iBinPsi) + 0.00001;
1933 binPsiUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->GetBinUpEdge(iBinPsi) - 0.00001;
1935 binVertexLowEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinLowEdge(iBinVertex) + 0.00001;
1936 binVertexUpEdge = fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->GetBinUpEdge(iBinVertex) - 0.00001;
1939 // get the 2D histograms for the correct type
1941 fSame = GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1942 fMixed = bMixed->GetCorrelationFunctionPN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1944 else if(type=="NP"){
1945 fSame = GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1946 fMixed = bMixed->GetCorrelationFunctionNP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1948 else if(type=="PP"){
1949 fSame = GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1950 fMixed = bMixed->GetCorrelationFunctionPP(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1952 else if(type=="NN"){
1953 fSame = GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1954 fMixed = bMixed->GetCorrelationFunctionNN(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1956 else if(type=="ALL"){
1957 fSame = GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1958 fMixed = bMixed->GetCorrelationFunctionChargeIndependent(binPsiLowEdge,binPsiUpEdge,binVertexLowEdge,binVertexUpEdge,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1961 if(fMixed && normToTrig && fMixed->Integral()>0){
1963 // normalization of Event mixing to 1 at (0,0) --> Jan Fietes method
1964 // do it only on away-side (due to two-track cuts): pi +- pi/6.
1965 Int_t binXmin = fMixed->GetXaxis()->FindBin(0-10e-5);
1966 Int_t binXmax = fMixed->GetXaxis()->FindBin(0+10e-5);
1967 Double_t binsX = (Double_t)(binXmax - binXmin + 1);
1968 Int_t binYmin = fMixed->GetYaxis()->FindBin(TMath::Pi() - normalizationRangePhi);
1969 Int_t binYmax = fMixed->GetYaxis()->FindBin(TMath::Pi() + normalizationRangePhi - 0.00001);
1970 Double_t binsY = (Double_t)(binYmax - binYmin + 1);
1972 Double_t mixedNorm = fMixed->Integral(binXmin,binXmax,binYmin,binYmax);
1973 mixedNorm /= binsX * binsY;
1975 // finite bin correction
1976 Double_t binWidthEta = fMixed->GetXaxis()->GetBinWidth(fMixed->GetNbinsX());
1977 Double_t maxEta = fMixed->GetXaxis()->GetBinUpEdge(fMixed->GetNbinsX());
1979 Double_t finiteBinCorrection = -1.0 / (2*maxEta) * binWidthEta / 2 + 1;
1980 //Printf("Finite bin correction: %f", finiteBinCorrection);
1981 mixedNorm /= finiteBinCorrection;
1983 fMixed->Scale(1./mixedNorm);
1986 if(fSame && fMixed){
1987 // then get the correlation function (divide fSame/fmixed)
1989 if(fMixed->Integral()>0)
1990 fSame->Divide(fMixed);
1992 // averaging with number of triggers:
1993 // average over number of triggers in each sub-bin
1994 Double_t NTrigSubBin = 0;
1995 if(type=="PN" || type=="PP")
1996 NTrigSubBin = (Double_t)(fHistP->Project(0,1)->Integral());
1997 else if(type=="NP" || type=="NN")
1998 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral());
1999 else if(type=="ALL")
2000 NTrigSubBin = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
2001 fSame->Scale(NTrigSubBin);
2003 // only if event mixing has enough statistics
2004 if(fMixed->Integral()>0){
2006 // for the first: clone
2007 if( (iBinPsi == binPsiMin && iBinVertex == binVertexMin) || !gHist ){
2008 gHist = (TH2D*)fSame->Clone();
2010 else{ // otherwise: add for averaging
2015 // otherwise record the number of triggers for correct normalization
2017 nTrigSubBinEmpty += NTrigSubBin;
2026 // averaging with number of triggers:
2027 // first set to full range and then obtain number of all triggers
2028 Double_t NTrigAll = 0;
2029 if(type=="PN" || type=="PP"){
2030 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2031 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2032 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2033 NTrigAll = (Double_t)(fHistP->Project(0,1)->Integral());
2035 else if(type=="NP" || type=="NN"){
2036 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2037 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2038 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2039 NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral());
2041 else if(type=="ALL"){
2042 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2043 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2044 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2045 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2046 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2047 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2048 NTrigAll = (Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral());
2051 // subtract number of triggers with empty sub bins for correct normalization
2052 NTrigAll -= nTrigSubBinEmpty;
2054 gHist->Scale(1./NTrigAll);
2061 //____________________________________________________________________//
2062 TH2D *AliBalancePsi::GetCorrelationFunctionPN(Double_t psiMin,
2064 Double_t vertexZMin,
2065 Double_t vertexZMax,
2066 Double_t ptTriggerMin,
2067 Double_t ptTriggerMax,
2068 Double_t ptAssociatedMin,
2069 Double_t ptAssociatedMax) {
2070 //Returns the 2D correlation function for (+-) pairs
2073 if(psiMin > psiMax-0.00001){
2074 AliError("psiMax <= psiMin");
2077 if(vertexZMin > vertexZMax-0.00001){
2078 AliError("vertexZMax <= vertexZMin");
2081 if(ptTriggerMin > ptTriggerMax-0.00001){
2082 AliError("ptTriggerMax <= ptTriggerMin");
2085 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2086 AliError("ptAssociatedMax <= ptAssociatedMin");
2091 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2092 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2096 //Printf("Cutting on Vz...");
2097 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2098 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2102 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2103 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2104 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2108 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2109 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2111 //fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2112 //fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(-0.5,2.5);
2114 //TH2D *gHistTest = dynamic_cast<TH2D *>(fHistP->Project(0,0,1));
2115 //TCanvas *c1 = new TCanvas("c1","");
2118 //AliError("Projection of fHistP = NULL");
2122 //gHistTest->DrawCopy("colz");
2125 //0:step, 1: Delta eta, 2: Delta phi
2126 TH2D *gHist = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2128 AliError("Projection of fHistPN = NULL");
2132 //AliInfo(Form("Entries (test): %lf",(Double_t)(gHistTest->GetEntries())));
2133 //AliInfo(Form("Entries (1D): %lf",(Double_t)(fHistP->Project(0,1)->GetEntries())));
2134 //AliInfo(Form("Entries (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->GetEntries())));
2135 //AliInfo(Form("Integral (1D): %lf",(Double_t)(fHistP->Project(0,1)->Integral())));
2136 //AliInfo(Form("Integral (2D): %lf",(Double_t)(fHistPN->Project(0,1,2)->Integral())));
2138 //TCanvas *c2 = new TCanvas("c2","");
2140 //fHistPN->Project(0,1,2)->DrawCopy("colz");
2142 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2143 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2145 //normalize to bin width
2146 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2151 //____________________________________________________________________//
2152 TH2D *AliBalancePsi::GetCorrelationFunctionNP(Double_t psiMin,
2154 Double_t vertexZMin,
2155 Double_t vertexZMax,
2156 Double_t ptTriggerMin,
2157 Double_t ptTriggerMax,
2158 Double_t ptAssociatedMin,
2159 Double_t ptAssociatedMax) {
2160 //Returns the 2D correlation function for (+-) pairs
2163 if(psiMin > psiMax-0.00001){
2164 AliError("psiMax <= psiMin");
2167 if(vertexZMin > vertexZMax-0.00001){
2168 AliError("vertexZMax <= vertexZMin");
2171 if(ptTriggerMin > ptTriggerMax-0.00001){
2172 AliError("ptTriggerMax <= ptTriggerMin");
2175 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2176 AliError("ptAssociatedMax <= ptAssociatedMin");
2181 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2182 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2186 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2187 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2191 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2192 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2193 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2197 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2198 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2200 //0:step, 1: Delta eta, 2: Delta phi
2201 TH2D *gHist = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2203 AliError("Projection of fHistPN = NULL");
2207 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2208 //Printf("Entries (2D): %lf",(Double_t)(fHistNP->Project(0,2,3)->GetEntries()));
2209 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2210 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2212 //normalize to bin width
2213 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2218 //____________________________________________________________________//
2219 TH2D *AliBalancePsi::GetCorrelationFunctionPP(Double_t psiMin,
2221 Double_t vertexZMin,
2222 Double_t vertexZMax,
2223 Double_t ptTriggerMin,
2224 Double_t ptTriggerMax,
2225 Double_t ptAssociatedMin,
2226 Double_t ptAssociatedMax) {
2227 //Returns the 2D correlation function for (+-) pairs
2230 if(psiMin > psiMax-0.00001){
2231 AliError("psiMax <= psiMin");
2234 if(vertexZMin > vertexZMax-0.00001){
2235 AliError("vertexZMax <= vertexZMin");
2238 if(ptTriggerMin > ptTriggerMax-0.00001){
2239 AliError("ptTriggerMax <= ptTriggerMin");
2242 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2243 AliError("ptAssociatedMax <= ptAssociatedMin");
2248 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2249 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2253 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2254 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2258 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2259 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2260 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2264 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2265 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2267 //0:step, 1: Delta eta, 2: Delta phi
2268 TH2D *gHist = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2270 AliError("Projection of fHistPN = NULL");
2274 //Printf("Entries (1D): %lf",(Double_t)(fHistP->Project(0,2)->GetEntries()));
2275 //Printf("Entries (2D): %lf",(Double_t)(fHistPP->Project(0,2,3)->GetEntries()));
2276 if((Double_t)(fHistP->Project(0,1)->Integral())>0)
2277 gHist->Scale(1./(Double_t)(fHistP->Project(0,1)->Integral()));
2279 //normalize to bin width
2280 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2285 //____________________________________________________________________//
2286 TH2D *AliBalancePsi::GetCorrelationFunctionNN(Double_t psiMin,
2288 Double_t vertexZMin,
2289 Double_t vertexZMax,
2290 Double_t ptTriggerMin,
2291 Double_t ptTriggerMax,
2292 Double_t ptAssociatedMin,
2293 Double_t ptAssociatedMax) {
2294 //Returns the 2D correlation function for (+-) pairs
2297 if(psiMin > psiMax-0.00001){
2298 AliError("psiMax <= psiMin");
2301 if(vertexZMin > vertexZMax-0.00001){
2302 AliError("vertexZMax <= vertexZMin");
2305 if(ptTriggerMin > ptTriggerMax-0.00001){
2306 AliError("ptTriggerMax <= ptTriggerMin");
2309 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2310 AliError("ptAssociatedMax <= ptAssociatedMin");
2315 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2316 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2320 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2321 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2325 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2326 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2327 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2331 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.))
2332 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2334 //0:step, 1: Delta eta, 2: Delta phi
2335 TH2D *gHist = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2337 AliError("Projection of fHistPN = NULL");
2341 //Printf("Entries (1D): %lf",(Double_t)(fHistN->Project(0,2)->GetEntries()));
2342 //Printf("Entries (2D): %lf",(Double_t)(fHistNN->Project(0,2,3)->GetEntries()));
2343 if((Double_t)(fHistN->Project(0,1)->Integral())>0)
2344 gHist->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral()));
2346 //normalize to bin width
2347 gHist->Scale(1./((Double_t)gHist->GetXaxis()->GetBinWidth(1)*(Double_t)gHist->GetYaxis()->GetBinWidth(1)));
2352 //____________________________________________________________________//
2353 TH2D *AliBalancePsi::GetCorrelationFunctionChargeIndependent(Double_t psiMin,
2355 Double_t vertexZMin,
2356 Double_t vertexZMax,
2357 Double_t ptTriggerMin,
2358 Double_t ptTriggerMax,
2359 Double_t ptAssociatedMin,
2360 Double_t ptAssociatedMax) {
2361 //Returns the 2D correlation function for the sum of all charge combination pairs
2364 if(psiMin > psiMax-0.00001){
2365 AliError("psiMax <= psiMin");
2368 if(vertexZMin > vertexZMax-0.00001){
2369 AliError("vertexZMax <= vertexZMin");
2372 if(ptTriggerMin > ptTriggerMax-0.00001){
2373 AliError("ptTriggerMax <= ptTriggerMin");
2376 if(ptAssociatedMin > ptAssociatedMax-0.00001){
2377 AliError("ptAssociatedMax <= ptAssociatedMin");
2382 fHistN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2383 fHistP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2384 fHistNN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2385 fHistPP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2386 fHistNP->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2387 fHistPN->GetGrid(0)->GetGrid()->GetAxis(0)->SetRangeUser(psiMin,psiMax-0.00001);
2391 fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2392 fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2393 fHistPN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2394 fHistNP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2395 fHistPP->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2396 fHistNN->GetGrid(0)->GetGrid()->GetAxis(5)->SetRangeUser(vertexZMin,vertexZMax-0.00001);
2400 if((ptTriggerMin != -1.)&&(ptTriggerMax != -1.)) {
2401 fHistN->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2402 fHistP->GetGrid(0)->GetGrid()->GetAxis(1)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2403 fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2404 fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2405 fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2406 fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRangeUser(ptTriggerMin,ptTriggerMax-0.00001);
2410 if((ptAssociatedMin != -1.)&&(ptAssociatedMax != -1.)){
2411 fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2412 fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2413 fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2414 fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRangeUser(ptAssociatedMin,ptAssociatedMax-0.00001);
2417 //0:step, 1: Delta eta, 2: Delta phi
2418 TH2D *gHistNN = dynamic_cast<TH2D *>(fHistNN->Project(0,1,2));
2420 AliError("Projection of fHistNN = NULL");
2423 TH2D *gHistPP = dynamic_cast<TH2D *>(fHistPP->Project(0,1,2));
2425 AliError("Projection of fHistPP = NULL");
2428 TH2D *gHistNP = dynamic_cast<TH2D *>(fHistNP->Project(0,1,2));
2430 AliError("Projection of fHistNP = NULL");
2433 TH2D *gHistPN = dynamic_cast<TH2D *>(fHistPN->Project(0,1,2));
2435 AliError("Projection of fHistPN = NULL");
2439 // sum all 2 particle histograms
2440 gHistNN->Add(gHistPP);
2441 gHistNN->Add(gHistNP);
2442 gHistNN->Add(gHistPN);
2444 // divide by sum of + and - triggers
2445 if((Double_t)(fHistN->Project(0,1)->Integral())>0 && (Double_t)(fHistP->Project(0,1)->Integral())>0)
2446 gHistNN->Scale(1./(Double_t)(fHistN->Project(0,1)->Integral() + fHistP->Project(0,1)->Integral()));
2448 //normalize to bin width
2449 gHistNN->Scale(1./((Double_t)gHistNN->GetXaxis()->GetBinWidth(1)*(Double_t)gHistNN->GetYaxis()->GetBinWidth(1)));
2454 //____________________________________________________________________//
2455 Bool_t AliBalancePsi::GetMomentsAnalytical(Int_t fVariable, TH1D* gHist, Bool_t kUseZYAM,
2456 Double_t &mean, Double_t &meanError,
2457 Double_t &sigma, Double_t &sigmaError,
2458 Double_t &skewness, Double_t &skewnessError,
2459 Double_t &kurtosis, Double_t &kurtosisError) {
2461 // helper method to calculate the moments and errors of a TH1D anlytically
2462 // fVariable = 1 for Delta eta (moments in full range)
2463 // = 2 for Delta phi (moments only on near side -pi/2 < dphi < pi/2)
2466 Bool_t success = kFALSE;
2472 skewnessError = -1.;
2474 kurtosisError = -1.;
2478 // ----------------------------------------------------------------------
2479 // basic parameters of histogram
2481 Int_t fNumberOfBins = gHist->GetNbinsX();
2482 // Int_t fBinWidth = gHist->GetBinWidth(1); // assume equal binning
2485 // ----------------------------------------------------------------------
2486 // ZYAM (for partially negative distributions)
2487 // --> we subtract always the minimum value
2488 Double_t zeroYield = 0.;
2489 Double_t zeroYieldCur = -FLT_MAX;
2491 for(Int_t iMin = 0; iMin<2; iMin++){
2492 zeroYieldCur = gHist->GetMinimum(zeroYieldCur);
2493 zeroYield += zeroYieldCur;
2496 //zeroYield = gHist->GetMinimum();
2498 // ----------------------------------------------------------------------
2499 // first calculate the mean
2501 Double_t fWeightedAverage = 0.;
2502 Double_t fNormalization = 0.;
2504 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2506 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2507 if(fVariable == 2 &&
2508 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2509 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2513 fWeightedAverage += (gHist->GetBinContent(i)-zeroYield) * gHist->GetBinCenter(i);
2514 fNormalization += (gHist->GetBinContent(i)-zeroYield);
2517 mean = fWeightedAverage / fNormalization;
2519 // ----------------------------------------------------------------------
2520 // then calculate the higher moments
2531 for(Int_t i = 1; i <= fNumberOfBins; i++) {
2533 // for Delta phi: moments only on near side -pi/2 < dphi < pi/2
2534 if(fVariable == 2 &&
2535 (gHist->GetBinCenter(i) < - TMath::Pi()/2 ||
2536 gHist->GetBinCenter(i) > TMath::Pi()/2)){
2540 fMu += (gHist->GetBinContent(i)-zeroYield) * (gHist->GetBinCenter(i) - mean);
2541 fMu2 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),2);
2542 fMu3 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),3);
2543 fMu4 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),4);
2544 fMu5 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),5);
2545 fMu6 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),6);
2546 fMu7 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),7);
2547 fMu8 += (gHist->GetBinContent(i)-zeroYield) * TMath::Power((gHist->GetBinCenter(i) - mean),8);
2550 // normalize to bin entries!
2551 fMu /= fNormalization;
2552 fMu2 /= fNormalization;
2553 fMu3 /= fNormalization;
2554 fMu4 /= fNormalization;
2555 fMu5 /= fNormalization;
2556 fMu6 /= fNormalization;
2557 fMu7 /= fNormalization;
2558 fMu8 /= fNormalization;
2560 sigma = TMath::Sqrt(fMu2);
2561 skewness = fMu3 / TMath::Power(sigma,3);
2562 kurtosis = fMu4 / TMath::Power(sigma,4) - 3;
2564 // normalize with sigma^r
2565 fMu /= TMath::Power(sigma,1);
2566 fMu2 /= TMath::Power(sigma,2);
2567 fMu3 /= TMath::Power(sigma,3);
2568 fMu4 /= TMath::Power(sigma,4);
2569 fMu5 /= TMath::Power(sigma,5);
2570 fMu6 /= TMath::Power(sigma,6);
2571 fMu7 /= TMath::Power(sigma,7);
2572 fMu8 /= TMath::Power(sigma,8);
2574 // ----------------------------------------------------------------------
2575 // then calculate the higher moment errors
2576 // cout<<fNormalization<<" "<<gHist->GetEffectiveEntries()<<" "<<gHist->Integral()<<endl;
2577 // cout<<gHist->GetMean()<<" "<<gHist->GetRMS()<<" "<<gHist->GetSkewness(1)<<" "<<gHist->GetKurtosis(1)<<endl;
2578 // cout<<gHist->GetMeanError()<<" "<<gHist->GetRMSError()<<" "<<gHist->GetSkewness(11)<<" "<<gHist->GetKurtosis(11)<<endl;
2580 Double_t normError = gHist->GetEffectiveEntries(); //gives the "ROOT" meanError
2582 // // assuming normal distribution (as done in ROOT)
2583 // meanError = sigma / TMath::Sqrt(normError);
2584 // sigmaError = TMath::Sqrt(sigma*sigma/(2*normError));
2585 // skewnessError = TMath::Sqrt(6./(normError));
2586 // kurtosisError = TMath::Sqrt(24./(normError));
2588 // use delta theorem paper (Luo - arXiv:1109.0593v1)
2589 Double_t Lambda11 = TMath::Abs((fMu4-1)*sigma*sigma/(4*normError));
2590 Double_t Lambda22 = TMath::Abs((9-6*fMu4+fMu3*fMu3*(35+9*fMu4)/4-3*fMu3*fMu5+fMu6)/normError);
2591 Double_t Lambda33 = TMath::Abs((-fMu4*fMu4+4*fMu4*fMu4*fMu4+16*fMu3*fMu3*(1+fMu4)-8*fMu3*fMu5-4*fMu4*fMu6+fMu8)/normError);
2592 //Double_t Lambda12 = -(fMu3*(5+3*fMu4)-2*fMu5)*sigma/(4*normError);
2593 //Double_t Lambda13 = ((-4*fMu3*fMu3+fMu4-2*fMu4*fMu4+fMu6)*sigma)/(2*normError);
2594 //Double_t Lambda23 = (6*fMu3*fMu3*fMu3-(3+2*fMu4)*fMu5+3*fMu3*(8+fMu4+2*fMu4*fMu4-fMu6)/2+fMu7)/normError;
2596 // cout<<Lambda11<<" "<<Lambda22<<" "<<Lambda33<<" "<<endl;
2597 // cout<<Lambda12<<" "<<Lambda13<<" "<<Lambda23<<" "<<endl;
2599 if (TMath::Sqrt(normError) != 0){
2600 meanError = sigma / TMath::Sqrt(normError);
2601 sigmaError = TMath::Sqrt(Lambda11);
2602 skewnessError = TMath::Sqrt(Lambda22);
2603 kurtosisError = TMath::Sqrt(Lambda33);
2608 else success = kFALSE;
2615 //____________________________________________________________________//
2616 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) {
2618 // calculates dphistar
2620 Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2622 static const Double_t kPi = TMath::Pi();
2625 // if (dphistar > 2 * kPi)
2626 // dphistar -= 2 * kPi;
2627 // if (dphistar < -2 * kPi)
2628 // dphistar += 2 * kPi;
2631 dphistar = kPi * 2 - dphistar;
2632 if (dphistar < -kPi)
2633 dphistar = -kPi * 2 - dphistar;
2634 if (dphistar > kPi) // might look funny but is needed
2635 dphistar = kPi * 2 - dphistar;
2640 //____________________________________________________________________//
2641 Double_t* AliBalancePsi::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2643 // This method is a copy from AliUEHist::GetBinning
2644 // takes the binning from <configuration> identified by <tag>
2645 // configuration syntax example:
2646 // 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
2649 // returns bin edges which have to be deleted by the caller
2651 TString config(configuration);
2652 TObjArray* lines = config.Tokenize("\n");
2653 for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2655 TString line(lines->At(i)->GetName());
2656 if (line.BeginsWith(TString(tag) + ":"))
2658 line.Remove(0, strlen(tag) + 1);
2659 line.ReplaceAll(" ", "");
2660 TObjArray* binning = line.Tokenize(",");
2661 Double_t* bins = new Double_t[binning->GetEntriesFast()];
2662 for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2663 bins[j] = TString(binning->At(j)->GetName()).Atof();
2665 nBins = binning->GetEntriesFast() - 1;
2674 AliFatal(Form("Tag %s not found in %s", tag, configuration));
2678 //____________________________________________________________________//
2679 Double_t AliBalancePsi::GetFWHM(Int_t gDeltaEtaPhi, TH1D* gHist,
2680 Double_t &fwhm_spline, Double_t &fwhmError) {
2682 // helper method to calculate the fwhm and errors of a TH1D
2684 Int_t repeat = 10000;
2686 if (gDeltaEtaPhi == 1){
2691 Double_t xmin = gHist->GetXaxis()->GetXmin();
2692 Double_t xmax = gHist->GetXaxis()->GetXmax();
2694 //+++++++FWHM TSpline+++++++++++++++++++++++++++//
2695 TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0);
2696 spline3->SetLineColor(kGreen+2);
2697 spline3->SetLineWidth(2);
2698 //spline3->Draw("SAME");
2700 Int_t Nbin = gHist->GetNbinsX();
2701 Int_t bin_max_hist_y = gHist->GetMaximumBin();
2702 Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
2703 Double_t max_spline = spline3->Eval(bins_center_x);
2704 Double_t y_spline = -999;
2705 Double_t y_spline_r = -999;
2706 Double_t x_spline_l = -999;
2707 Double_t x_spline_r = -999;
2709 for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
2710 y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
2711 if (y_spline > max_spline/2){
2712 x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
2716 for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
2717 y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
2718 if (y_spline_r < max_spline/2){
2719 x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
2723 fwhm_spline = x_spline_r - x_spline_l;
2724 //+++++++FWHM TSpline++++++++++++++++++++++++//
2726 //+++++++Error Calculation SPLINE++++++++++++//
2727 Double_t dmeans,dsigmas;
2729 Int_t bin_max_hist_y1;
2730 Double_t bins_center_x1;
2731 Double_t max_spline1;
2732 Double_t y_spline1 = - 999.;
2733 Double_t y_spline_r1 = - 999.;
2734 Double_t x_spline_l1 = -999.;
2735 Double_t x_spline_r1 = -999.;
2736 Double_t fwhm_spline1 = 0.;
2737 Double_t fwhm_T = 0.;
2739 TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
2740 TRandom3 *rgauss = new TRandom3(0);
2741 TH1F *hists = new TH1F("hists","hists",1000,1.,3);
2743 for (Int_t f=0; f<repeat; f++){ //100
2744 for (Int_t h=0; h<Nbin+1; h++){
2745 dmeans = gHist->GetBinContent(h);
2746 dsigmas = gHist->GetBinError(h);
2747 Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
2748 hEmpty2->SetBinContent(h,rgauss_point);
2754 spline1 = new TSpline3(hEmpty2,"b2e2",0,0);
2755 spline1->SetLineColor(kMagenta+1);
2756 spline1->SetLineWidth(1);
2757 //spline1->Draw("SAME");
2759 bin_max_hist_y1 = hEmpty2->GetMaximumBin();
2760 bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
2761 max_spline1 = spline1->Eval(bins_center_x1);
2763 for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
2764 y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
2765 if (y_spline1 > max_spline1/2){
2766 x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
2770 for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
2771 y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
2772 if (y_spline_r1 < max_spline1/2){
2773 x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
2778 fwhm_spline1 = x_spline_r1 - x_spline_l1;
2779 hists->Fill(fwhm_spline1);
2780 fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);
2783 fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));
2784 //+++++++Error Calculation SPLINE+++++++++++//
2791 Double_t xmin = gHist->GetXaxis()->GetXmin();
2792 Double_t xmax = gHist->GetXaxis()->GetXmax();
2794 //+++++++FWHM TSpline+++++++++++++++++++++++++++//
2795 TSpline3 *spline3 = new TSpline3(gHist,"b2e2",0,0);
2796 spline3->SetLineColor(kGreen+2);
2797 spline3->SetLineWidth(2);
2798 //spline3->Draw("SAME");
2800 Int_t Nbin = gHist->GetNbinsX();
2801 Int_t bin_max_hist_y = gHist->GetMaximumBin();
2802 Double_t bins_center_x = gHist->GetBinCenter(bin_max_hist_y);
2803 Double_t max_spline = spline3->Eval(bins_center_x);
2805 Int_t bin_min_hist_y = gHist->GetMinimumBin();
2806 Double_t bins_center_xmin = gHist->GetBinCenter(bin_min_hist_y);
2807 Double_t min_spline = spline3->Eval(bins_center_xmin);
2809 Double_t y_spline = -999.;
2810 Double_t y_spline_r = -999.;
2811 Double_t x_spline_l = -999.;
2812 Double_t x_spline_r = -999.;
2814 for (Int_t i= 0; i < bin_max_hist_y*1000; i++){
2815 y_spline = spline3->Eval(xmin + i*(bins_center_x - xmin)/ (bin_max_hist_y*1000));
2816 if (y_spline > (max_spline+min_spline)/2){
2817 x_spline_l = xmin + (i-1)*(bins_center_x - xmin)/(bin_max_hist_y*1000);
2821 for (Int_t j= 0; j < bin_max_hist_y*1000; j++){
2822 y_spline_r = spline3->Eval(bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000));
2823 if (y_spline_r < (max_spline+min_spline)/2){
2824 x_spline_r = bins_center_x + j*(xmax - bins_center_x)/(bin_max_hist_y*1000);
2828 fwhm_spline = x_spline_r - x_spline_l;
2829 //+++++++FWHM TSpline++++++++++++++++++++++++//
2831 //+++++++Error Calculation SPLINE++++++++++++//
2832 Double_t dmeans,dsigmas;
2834 Int_t bin_max_hist_y1;
2835 Double_t bins_center_x1;
2836 Double_t max_spline1;
2837 Double_t y_spline1 = -999.;
2838 Double_t y_spline_r1 = -999.;
2839 Double_t x_spline_l1 = -999.;
2840 Double_t x_spline_r1 = -999.;
2841 Double_t fwhm_spline1 = 0.;
2842 Double_t fwhm_T = 0.;
2844 TH1F *hEmpty2 = new TH1F("hEmpty2","hEmpty2",Nbin,xmin,xmax);
2845 TRandom3 *rgauss = new TRandom3(0);
2846 Int_t bin_min_hist_y1;
2847 Double_t bins_center_xmin1,min_spline1;
2849 for (Int_t f=0; f<repeat; f++){ //100
2850 for (Int_t h=0; h<Nbin+1; h++){
2851 dmeans = gHist->GetBinContent(h);
2852 dsigmas = gHist->GetBinError(h);
2853 Double_t rgauss_point = rgauss->Gaus(dmeans,dsigmas);
2854 hEmpty2->SetBinContent(h,rgauss_point);
2860 spline1 = new TSpline3(hEmpty2,"b2e2",0,0);
2861 spline1->SetLineColor(kMagenta+1);
2862 spline1->SetLineWidth(1);
2863 //spline1->Draw("SAME");
2865 bin_max_hist_y1 = hEmpty2->GetMaximumBin();
2866 bins_center_x1 = hEmpty2->GetBinCenter(bin_max_hist_y1);
2867 max_spline1 = spline1->Eval(bins_center_x1);
2869 bin_min_hist_y1 = hEmpty2->GetMinimumBin();
2870 bins_center_xmin1 = hEmpty2->GetBinCenter(bin_min_hist_y1);
2871 min_spline1 = spline1->Eval(bins_center_xmin1);
2873 for (Int_t i= 0; i < bin_max_hist_y1*1000; i++){
2874 y_spline1 = spline3->Eval(xmin + i*(bins_center_x1 - xmin)/ (bin_max_hist_y1*1000));
2875 if (y_spline1 > (max_spline1+min_spline1)/2){
2876 x_spline_l1 = xmin + (i-1)*(bins_center_x1 - xmin)/(bin_max_hist_y1*1000);
2880 for (Int_t j= 0; j < bin_max_hist_y1*1000; j++){
2881 y_spline_r1 = spline3->Eval(bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000));
2882 if (y_spline_r1 < (max_spline1+min_spline1)/2){
2883 x_spline_r1 = bins_center_x1 + j*(xmax - bins_center_x1)/(bin_max_hist_y1*1000);
2888 fwhm_spline1 = x_spline_r1 - x_spline_l1;
2889 fwhm_T += (fwhm_spline - fwhm_spline1)*(fwhm_spline - fwhm_spline1);
2891 fwhmError = TMath::Sqrt(TMath::Abs(fwhm_T/(repeat-1)));