Event Mixing for Triggered BF analysis
[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalanceTriggered.cxx
1 /**************************************************************************
2  * Author: Panos Christakoglou.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 //-----------------------------------------------------------------
15 //           Balance Function class
16 //   This is the class to deal with the Balance Function analysis
17 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
18 //   Modified: Michael Weber, m.weber@cern.ch
19 //-----------------------------------------------------------------
20
21
22 //ROOT
23 #include <Riostream.h>
24 #include <TMath.h>
25 #include <TAxis.h>
26 #include <TH1D.h>
27 #include <TH2D.h>
28
29
30 #include <AliTHn.h>
31 #include <TLorentzVector.h>
32 #include <TObjArray.h>
33 #include <TGraphErrors.h>
34 #include <TString.h>
35
36 #include "AliVParticle.h"
37 #include "AliMCParticle.h"
38 #include "AliESDtrack.h"
39 #include "AliAODTrack.h"
40
41 #include "AliBalanceTriggered.h"
42
43 ClassImp(AliBalanceTriggered)
44
45 //____________________________________________________________________//
46 AliBalanceTriggered::AliBalanceTriggered() :
47   TObject(), 
48   fAnalysisLevel("AOD"),
49   fHistP(0x0),
50   fHistN(0x0),
51   fHistPN(0x0),
52   fHistNP(0x0),
53   fHistPP(0x0),
54   fHistNN(0x0)
55 {
56   // Default constructor
57  
58 }
59
60
61 //____________________________________________________________________//
62 AliBalanceTriggered::AliBalanceTriggered(const AliBalanceTriggered& balance):
63   TObject(balance), 
64   fAnalysisLevel(balance.fAnalysisLevel),
65   fHistP(balance.fHistP),
66   fHistN(balance.fHistN),
67   fHistPN(balance.fHistPN),
68   fHistNP(balance.fHistNP),
69   fHistPP(balance.fHistPP),
70   fHistNN(balance.fHistNN)
71 {
72   //copy constructor
73
74 }
75
76 //____________________________________________________________________//
77 AliBalanceTriggered::~AliBalanceTriggered() {
78   // Destructor
79  
80   delete fHistP;
81   delete fHistN;
82   delete fHistPN;
83   delete fHistNP;
84   delete fHistPP;
85   delete fHistNN;
86   
87
88 }
89
90 //____________________________________________________________________//
91 void AliBalanceTriggered::InitHistograms() {
92
93   //Initialize the histograms
94
95   TString title    = "";      // histogram title
96   Int_t anaSteps   = 1;       // analysis steps
97
98   // single particle histograms
99   Int_t iBinSingle[nTrackVarsSingle];         // binning for track variables
100   Double_t* dBinsSingle[nTrackVarsSingle];    // bins for track variables  
101   TString axisTitleSingle[nTrackVarsSingle];  // axis titles for track variables
102   
103   // two particle histograms
104   Int_t iBinPair[nTrackVarsPair];         // binning for track variables
105   Double_t* dBinsPair[nTrackVarsPair];    // bins for track variables  
106   TString axisTitlePair[nTrackVarsPair];  // axis titles for track variables
107
108    
109   //-----------------------------------------------------------
110   // histogram settings (hard coded!)
111   //-----------------------------------------------------------
112
113   // eta
114   const Int_t kNEtaBins = 40;
115   Double_t etaBins[kNEtaBins+1];
116   for(Int_t i = 0; i < kNEtaBins+1; i++){
117     etaBins[i] = -1.0 + i * 0.05;
118   } 
119   iBinSingle[0]       = kNEtaBins;
120   dBinsSingle[0]      = etaBins;
121   axisTitleSingle[0]  = "#eta"; 
122   
123   // delta eta
124   const Int_t kNDeltaEtaBins = 40;
125   Double_t deltaEtaBins[kNDeltaEtaBins+1];
126   for(Int_t i = 0; i < kNDeltaEtaBins+1; i++){
127     deltaEtaBins[i] = -2.0 + i * 0.1;
128   } 
129   iBinPair[0]       = kNDeltaEtaBins;
130   dBinsPair[0]      = deltaEtaBins;
131   axisTitlePair[0]  = "#Delta #eta"; 
132
133   // phi
134   const Int_t kNPhiBins = 72;
135   Double_t phiBins[kNPhiBins+1];
136   for(Int_t i = 0; i < kNPhiBins+1; i++){
137     phiBins[i] = 0.0 + i * 5.;
138   } 
139   iBinSingle[1]       = kNPhiBins;
140   dBinsSingle[1]      = phiBins;
141   axisTitleSingle[1]  = "#phi (#circ)"; 
142   
143   // delta phi
144   const Int_t kNDeltaPhiBins = 72;
145   Double_t deltaPhiBins[kNDeltaPhiBins+1];
146   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
147     deltaPhiBins[i] = -180.0 + i * 5.;
148   } 
149   iBinPair[1]       = kNDeltaPhiBins;
150   dBinsPair[1]      = deltaPhiBins;
151   axisTitlePair[1]  = "#Delta #phi (#circ)"; 
152
153   // pT
154   const Int_t kNPtBins = 22;
155   Double_t pTBins[kNPtBins+1] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
156   iBinSingle[2]       = kNPtBins;
157   dBinsSingle[2]      = pTBins;
158   axisTitleSingle[2]  = "p_{T,trig} (GeV/c)";
159  
160   iBinPair[2]       = kNPtBins;
161   dBinsPair[2]      = pTBins;
162   axisTitlePair[2]  = "p_{T} (GeV/c)";
163  
164   iBinPair[3]       = kNPtBins;
165   dBinsPair[3]      = pTBins;
166   axisTitlePair[3]  = "p_{T,trig} (GeV/c)"; 
167
168   // centrality
169   const Int_t kNCentBins = 9;
170   Double_t centBins[kNCentBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
171
172   iBinSingle[3]       = kNCentBins;
173   dBinsSingle[3]      = centBins;
174   axisTitleSingle[3]  = "centrality"; 
175
176   iBinPair[4]       = kNCentBins;
177   dBinsPair[4]      = centBins;
178   axisTitlePair[4]  = "centrality";
179
180   //-----------------------------------------------------------
181   
182
183
184   // User settings (depending on the analysis settings [only one at the moment])
185   if(fAnalysisLevel == "AOD"){
186     title = "hdEtaVsdPhi";
187   }
188
189   //-----------------------------------------------------------
190   //-----------------------------------------------------------
191   // Histogram creation
192
193   // histogram for negative particles
194   fHistN = new AliTHn(Form("fHistN"), Form("%s_N",title.Data()), anaSteps, nTrackVarsSingle, iBinSingle);
195   for (Int_t j=0; j<nTrackVarsSingle; j++)
196     {
197       fHistN->SetBinLimits(j, dBinsSingle[j]);
198       fHistN->SetVarTitle(j, axisTitleSingle[j]);
199     }
200
201   // histogram for positive particles
202   fHistP = new AliTHn(Form("fHistP"), Form("%s_P",title.Data()), anaSteps, nTrackVarsSingle, iBinSingle);
203   for (Int_t j=0; j<nTrackVarsSingle; j++)
204     {
205       fHistP->SetBinLimits(j, dBinsSingle[j]);
206       fHistP->SetVarTitle(j, axisTitleSingle[j]);
207     }
208
209   // histogram for +- pairs
210   fHistPN = new AliTHn(Form("fHistPN"), Form("%s_PN",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
211   for (Int_t j=0; j<nTrackVarsPair; j++)
212     {
213       fHistPN->SetBinLimits(j, dBinsPair[j]);
214       fHistPN->SetVarTitle(j, axisTitlePair[j]);
215     }
216
217   // histogram for -+ pairs
218   fHistNP = new AliTHn(Form("fHistNP"), Form("%s_NP",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
219   for (Int_t j=0; j<nTrackVarsPair; j++)
220     {
221       fHistNP->SetBinLimits(j, dBinsPair[j]);
222       fHistNP->SetVarTitle(j, axisTitlePair[j]);
223     }
224
225   // histogram for ++ pairs
226   fHistPP = new AliTHn(Form("fHistPP"), Form("%s_PP",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
227   for (Int_t j=0; j<nTrackVarsPair; j++)
228     {
229       fHistPP->SetBinLimits(j, dBinsPair[j]);
230       fHistPP->SetVarTitle(j, axisTitlePair[j]);
231     }
232
233   // histogram for -- pairs
234   fHistNN = new AliTHn(Form("fHistNN"), Form("%s_NN",title.Data()), anaSteps, nTrackVarsPair, iBinPair);
235   for (Int_t j=0; j<nTrackVarsPair; j++)
236     {
237       fHistNN->SetBinLimits(j, dBinsPair[j]);
238       fHistNN->SetVarTitle(j, axisTitlePair[j]);
239     }
240
241   //-----------------------------------------------------------
242   //-----------------------------------------------------------
243
244
245 }
246
247 //____________________________________________________________________//
248 void AliBalanceTriggered::FillBalance(Float_t fCentrality,TObjArray *particles, TObjArray *particlesMixed) {
249   // Calculates the balance function
250
251  
252   // Initialize histograms if not done yet
253   if(!fHistPN){
254     AliWarning("Histograms not yet initialized! --> Will be done now");
255     AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
256     InitHistograms();
257   }
258
259   Double_t trackVarsSingle[nTrackVarsSingle];
260   Double_t trackVarsPair[nTrackVarsPair];
261
262   if (!particles){
263     AliWarning("particles TObjArray is NULL pointer --> return");
264     return;
265   }
266   
267   // define end of particle loops
268   Int_t iMax = particles->GetEntriesFast();
269   Int_t jMax = iMax;
270   if (particlesMixed)
271     jMax = particlesMixed->GetEntriesFast();
272
273   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
274   TObjArray* particlesSecond = (particlesMixed) ? particlesMixed : particles;
275
276   TArrayF secondEta(jMax);
277   TArrayF secondPhi(jMax);
278   TArrayF secondPt(jMax);
279
280   for (Int_t i=0; i<jMax; i++){
281     secondEta[i] = ((AliVParticle*) particlesSecond->At(i))->Eta();
282     secondPhi[i] = ((AliVParticle*) particlesSecond->At(i))->Phi();
283     secondPt[i]  = ((AliVParticle*) particlesSecond->At(i))->Pt();
284   }
285   
286   // 1st particle loop
287   for (Int_t i=0; i<iMax; i++)
288     {
289       
290       AliVParticle* firstParticle = (AliVParticle*) particles->At(i);
291
292       // some optimization
293       Float_t firstEta = firstParticle->Eta();
294       Float_t firstPhi = firstParticle->Phi();
295       Float_t firstPt  = firstParticle->Pt();
296
297       
298       Short_t  charge = (Short_t) firstParticle->Charge();
299       trackVarsSingle[0]    =  firstEta;    //eta
300       trackVarsSingle[1]    =  firstPhi;    //phi
301       trackVarsSingle[2]    =  firstPt;     //pt trigger
302       trackVarsSingle[3]    =  fCentrality; //centrality 
303       
304       //fill single particle histograms
305       if(charge > 0)  fHistP->Fill(trackVarsSingle,0,1.); 
306       else            fHistN->Fill(trackVarsSingle,0,1.);  
307
308       
309       // 2nd particle loop (only for j < i for non double counting in the same pT region)
310       // --> SAME pT region for trigger and assoc: NO double counting with this
311       // --> DIFF pT region for trigger and assoc: Missing assoc. particles with j > i to a trigger i 
312       //                          --> can be handled afterwards by using assoc. as trigger as well ?!     
313       for(Int_t j = 0; j < i; j++) {   // or go to full here (everything prepared)?
314         
315         if (particlesMixed && jMax < i)  // if the mixed track number is smaller than the main event one (could be done better if one loops over all tracks)
316           break;
317
318         AliVParticle* secondParticle = (AliVParticle*) particlesSecond->At(j);
319
320         Short_t charge2 = (Short_t) secondParticle->Charge();
321         trackVarsPair[0]    =  firstEta - secondEta[j];  // delta eta
322         trackVarsPair[1]    =  firstPhi = secondPhi[j];  // delta phi
323         if (trackVarsPair[1] > 180)   // delta phi between -180 and 180 
324           trackVarsPair[1] -= 360;
325         if (trackVarsPair[1] <  - 180) 
326           trackVarsPair[1] += 360;
327         
328         trackVarsPair[2]    =  firstPt;      // pt trigger
329         trackVarsPair[3]    =  secondPt[j];  // pt
330         trackVarsPair[4]    =  fCentrality;  // centrality
331         
332         if( charge > 0 && charge2 < 0)  fHistPN->Fill(trackVarsPair,0,1.); 
333         else if( charge < 0 && charge2 > 0)  fHistNP->Fill(trackVarsPair,0,1.); 
334         else if( charge > 0 && charge2 > 0)  fHistPP->Fill(trackVarsPair,0,1.); 
335         else if( charge < 0 && charge2 < 0)  fHistNN->Fill(trackVarsPair,0,1.); 
336         else AliWarning("Wrong charge combination!");
337         
338       }//end of 2nd particle loop
339     }//end of 1st particle loop
340 }  
341
342
343 TH1D* AliBalanceTriggered::GetBalanceFunctionHistogram1D(Int_t var, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
344
345   // check which variable should be analyzed
346   // 0 = Delta eta
347   // 1 = Delta phi
348
349   if( var < 0 || var > 1){
350     AliError("Only Variable 0 (= Delta eta) or 1 (= Delta phi) allowed");
351     return NULL;
352   }
353
354
355   // Choose region to analyze 
356   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
357   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
358
359   // pT trigger
360   fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
361   fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
362   fHistPN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
363   fHistNP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
364   fHistPP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
365   fHistNN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
366
367   // pT
368   fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
369   fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
370   fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
371   fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
372
373   // centrality
374   fHistP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
375   fHistN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
376   fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
377   fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
378   fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
379   fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
380   
381
382   // Project into the wanted space (1st: analysis step, 2nd: axis)
383   TH1D* hTemp1 = (TH1D*)fHistPN->Project(0,var);
384   TH1D* hTemp2 = (TH1D*)fHistNP->Project(0,var);
385   TH1D* hTemp3 = (TH1D*)fHistPP->Project(0,var);
386   TH1D* hTemp4 = (TH1D*)fHistNN->Project(0,var);
387   TH1D* hTemp5 = (TH1D*)fHistP->Project(0,var);
388   TH1D* hTemp6 = (TH1D*)fHistN->Project(0,var);
389
390   TH1D* gHistBalanceFunctionHistogram = NULL;
391
392   // Calculate BF
393   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
394     
395     gHistBalanceFunctionHistogram = (TH1D*)hTemp1->Clone();
396     gHistBalanceFunctionHistogram->Reset();
397
398     // Calculate BF
399     hTemp1->Add(hTemp3,-1.);
400     hTemp1->Scale(1./hTemp5->GetEntries());
401     hTemp2->Add(hTemp4,-1.);
402     hTemp2->Scale(1./hTemp6->GetEntries());
403     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
404     gHistBalanceFunctionHistogram->Scale(0.5/1.);
405   }
406
407   return gHistBalanceFunctionHistogram;
408 }
409
410 TH2D* AliBalanceTriggered::GetBalanceFunctionHistogram2D(Int_t var1, Int_t var2, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
411
412   // check which variable should be analyzed
413   // 0 = Delta eta
414   // 1 = Delta phi
415
416   if( var1 < 0 || var1 > 1 || var2 < 0 || var2 > 1){
417     AliError("Only Variable 0 (= Delta eta) or 1 (= Delta phi) allowed");
418     return NULL;
419   }
420
421
422   // Choose region to analyze 
423   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
424   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
425
426   // pT trigger
427   fHistP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
428   fHistN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
429   fHistPN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
430   fHistNP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
431   fHistPP->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
432   fHistNN->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
433
434   // pT
435   fHistPN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
436   fHistNP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
437   fHistPP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
438   fHistNN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
439
440   // centrality
441   fHistP->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
442   fHistN->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
443   fHistPN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
444   fHistNP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
445   fHistPP->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
446   fHistNN->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
447   
448
449   // Project into the wanted space (1st: analysis step, 2nd: axis1, 3rd: axis2)
450   TH2D* hTemp1 = (TH2D*)fHistPN->Project(0,var1,var2);
451   TH2D* hTemp2 = (TH2D*)fHistNP->Project(0,var1,var2);
452   TH2D* hTemp3 = (TH2D*)fHistPP->Project(0,var1,var2);
453   TH2D* hTemp4 = (TH2D*)fHistNN->Project(0,var1,var2);
454   TH2D* hTemp5 = (TH2D*)fHistP->Project(0,var1,var2);
455   TH2D* hTemp6 = (TH2D*)fHistN->Project(0,var1,var2);
456
457   TH2D* gHistBalanceFunctionHistogram = NULL;
458
459   // Calculate BF
460   if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)&&(hTemp5)&&(hTemp6)) {
461     
462     gHistBalanceFunctionHistogram = (TH2D*)hTemp1->Clone();
463     gHistBalanceFunctionHistogram->Reset();
464     
465     // Calculate BF
466     hTemp1->Add(hTemp3,-1.);
467     hTemp1->Scale(1./hTemp5->GetEntries());
468     hTemp2->Add(hTemp4,-1.);
469     hTemp2->Scale(1./hTemp6->GetEntries());
470     gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
471   }
472
473   return gHistBalanceFunctionHistogram;
474 }
475
476
477 TH1D* AliBalanceTriggered::GetHistogram1D(Int_t histo, Int_t var, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
478
479   // check which variable should be analyzed
480   //
481   // pair histograms:
482   // 0 = Delta eta 
483   // 1 = Delta phi
484   // 2 = pT, trigger
485   // 3 = centrality
486   //
487   // pair histograms:
488   // 0 = Delta eta 
489   // 1 = Delta phi
490   // 2 = pT
491   // 3 = pT, trigger
492   // 4 = centrality
493
494   if(histo < 0 || histo > 5){
495     AliError("Only 6 histograms available: 0(P), 1(N), 2(PN), 3(NP), 4(PP), 5(NN)");
496     return NULL;
497   }
498
499   if( histo > 1 && (var < 0 || var > 5)){
500     AliError("Only Variable 0 to 4 allowed for pair histograms (histo > 1)");
501     return NULL;
502   }
503   if( histo < 2 && (var < 0 || var > 4)){
504     AliError("Only Variable 0 to 3 allowed for single histograms (histo < 2)");
505     return NULL;
506   }
507
508   // get the histogram
509   AliTHn *gTHn = NULL;
510   switch(histo){
511  
512   case 0:
513     gTHn = fHistP;
514     break;
515
516   case 1:
517     gTHn = fHistN;
518     break;
519
520   case 2:
521     gTHn = fHistPN;
522     break;
523
524   case 3:
525     gTHn = fHistNP;
526     break;
527
528   case 4:
529     gTHn = fHistPP;
530     break;
531
532   case 5:
533     gTHn = fHistNN;
534     break;
535     
536   }
537
538   if(!gTHn){
539     AliError(Form("AliTHn number %d = NULL",histo));
540     return NULL;
541   }
542
543   // Choose region to analyze 
544   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
545   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
546
547   // pT trigger
548   gTHn->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
549  
550   // pT
551   if(histo > 1) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
552  
553   // centrality
554   if(histo < 2) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
555   else          gTHn->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
556  
557
558   // Project into the wanted space (1st: analysis step, 2nd: axis)
559   TH1D* gHisto = (TH1D*)gTHn->Project(0,var);
560
561   return gHisto;
562 }
563
564
565 TH2D* AliBalanceTriggered::GetHistogram2D(Int_t histo, Int_t var1, Int_t var2, Double_t pTMinTrigger, Double_t pTMaxTrigger, Double_t pTMin, Double_t pTMax, Double_t centrMin, Double_t centrMax){
566
567   // check which variable should be analyzed
568   //
569   // pair histograms:
570   // 0 = Delta eta 
571   // 1 = Delta phi
572   // 2 = pT, trigger
573   // 3 = centrality
574   //
575   // pair histograms:
576   // 0 = Delta eta 
577   // 1 = Delta phi
578   // 2 = pT
579   // 3 = pT, trigger
580   // 4 = centrality
581
582   if(histo < 0 || histo > 5){
583     AliError("Only 6 histograms available: 0(P), 1(N), 2(PN), 3(NP), 4(PP), 5(NN)");
584     return NULL;
585   }
586
587   if( histo > 1 && (var1 < 0 || var1 > 5 || var2 < 0 || var2 > 5)){
588     AliError("Only Variable 0 to 4 allowed for pair histograms (histo > 1)");
589     return NULL;
590   }
591   if( histo < 2 && (var1 < 0 || var1 > 4 || var2 < 0 || var2 > 4)){
592     AliError("Only Variable 0 to 3 allowed for single histograms (histo < 2)");
593     return NULL;
594   }
595
596   // get the histogram
597   AliTHn *gTHn = NULL;
598   switch(histo){
599  
600   case 0:
601     gTHn = fHistP;
602     break;
603
604   case 1:
605     gTHn = fHistN;
606     break;
607
608   case 2:
609     gTHn = fHistPN;
610     break;
611
612   case 3:
613     gTHn = fHistNP;
614     break;
615
616   case 4:
617     gTHn = fHistPP;
618     break;
619
620   case 5:
621     gTHn = fHistNN;
622     break;
623     
624   }
625
626   if(!gTHn){
627     AliError(Form("AliTHn number %d = NULL",histo));
628     return NULL;
629   }
630
631   // Choose region to analyze 
632   // for Single Histograms (P,N):       2 = pT,trigger; 3 = centrality
633   // for Pair Histograms (PN,NP,NN,PP): 2 = pT; 3 = pT,trigger; 4 = centrality
634
635   // pT trigger
636   gTHn->GetGrid(0)->GetGrid()->GetAxis(2)->SetRange(pTMinTrigger,pTMaxTrigger); 
637  
638   // pT
639   if(histo > 1) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(pTMin,pTMax); 
640  
641   // centrality
642   if(histo < 2) gTHn->GetGrid(0)->GetGrid()->GetAxis(3)->SetRange(centrMin,centrMax); 
643   else          gTHn->GetGrid(0)->GetGrid()->GetAxis(4)->SetRange(centrMin,centrMax); 
644  
645
646   // Project into the wanted space (1st: analysis step, 2nd: axis)
647   TH2D* gHisto = (TH2D*)gTHn->Project(0,var1,var2);
648
649   return gHisto;
650 }