AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetJTJT.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 //
16 //
17 // Jet fragmentation transverse momentum (j_T) analysis task
18 //
19 // Author: T.Snellman
20 
21 #include <TClonesArray.h>
22 #include <TH1F.h>
23 #include <TH2F.h>
24 #include <TH3F.h>
25 #include <TList.h>
26 #include <TProfile.h>
27 #include <TLorentzVector.h>
28 #include <TVector.h>
29 #include <TGraphErrors.h>
30 #include <TGrid.h>
31 #include <TSystem.h>
32 #include <TFile.h>
33 
34 #include "AliCentrality.h"
35 
36 
37 
38 #include "AliVCluster.h"
39 #include "AliAODCaloCluster.h"
40 #include "AliESDCaloCluster.h"
41 #include "AliVTrack.h"
42 #include "AliAnalysisUtils.h"
43 #include "AliEmcalJet.h"
44 #include "AliRhoParameter.h"
45 #include "AliLog.h"
46 #include "AliJetContainer.h"
47 #include "AliParticleContainer.h"
48 #include "AliClusterContainer.h"
49 #include "AliPicoTrack.h"
50 
51 #include "AliAnalysisTaskJetJTJT.h"
52 
53 
54 ClassImp(AliAnalysisTaskJetJTJT)
55 
56  //________________________________________________________________________
59  fHistTracksPt(0),
60  fHistTracksJt(0),
61  fHistTracksEta(0),
62  fHistClustersPt(0),
63  fHistLeadingJetPt(0),
64  fHistJetsPt(0),
65  fHistJetsCorrPt(0),
66  fHistJetsCorrPtVsNonCorr(0),
67  fHistBackgroundDone(0),
68  fHistJTPta(0),
69  fHistLogJTPta(0),
70  fHistJTPta_all(0),
71  fHistJTBg(0),
72  fHistLogJTBg(0),
73  fHistPtaVsJt(0),
74  fHistBgPtaVsJt(0),
75  fHistJTPtaNonInv(0),
76  fHistLogJTPtaNonInv(0),
77  fHistJTPta_allNonInv(0),
78  fHistJTBgNonInv(0),
79  fHistLogJTBgNonInv(0),
80  fHistBgMulti(0),
81  fHistBgPt(0),
82  fHistJetEta(0),
83  fHistJetMulti(0),
84  fHistJetTracksPt(0),
85  fhTrackingEfficiency(0),
86  fNpttBins(1),
87  fNptaBins(1),
88  fEffMode(1),
89  fJetsCont(0),
90  //fJetsConts(0),
91  //nJetsConts(0),
92  fTracksCont(0),
93  fCaloClustersCont(0),
94  fhVertexZ(0),
95  fHistEvtSelection(0),
96  fPrimaryVertex(0),
97  fTracks(0),
98  fTrackArrayName("nonejk"),
99  runPeriod(""),
100  fEfficiency(0),
101  fVertexHelper(0),
102  debug(0)
103 
104 
105 {
106  // Default constructor.
107 
108  fHistTracksPt = new TH1*[fNcentBins];
109  fHistTracksJt = new TH1*[fNcentBins];
110  fHistTracksEta = new TH1*[fNcentBins];
111  fHistClustersPt = new TH1*[fNcentBins];
112  fHistLeadingJetPt = new TH1*[fNcentBins];
113  fHistJetsPt = new TH1**[fNcentBins];
114  fHistJetsCorrPt = new TH1**[fNcentBins];
115  fHistJetsCorrPtVsNonCorr= new TProfile*[fNcentBins];
116  fHistBackgroundDone = new TH1**[fNcentBins];
117  fHistJTPta = new TH1***[fNcentBins];
118  fHistLogJTPta = new TH1***[fNcentBins];
119  fHistJTPta_all = new TH1***[fNcentBins];
120  fHistJTBg = new TH1***[fNcentBins];
121  fHistLogJTBg = new TH1***[fNcentBins];
122  fHistJTPtaNonInv = new TH1***[fNcentBins];
123  fHistLogJTPtaNonInv = new TH1***[fNcentBins];
124  fHistJTPta_allNonInv = new TH1***[fNcentBins];
125  fHistJTBgNonInv = new TH1***[fNcentBins];
126  fHistLogJTBgNonInv = new TH1***[fNcentBins];
127  fHistBgMulti = new TH1**[fNcentBins];
128  fHistBgPt = new TH1**[fNcentBins];
129  fHistJetEta = new TH1**[fNcentBins];
130  fHistJetMulti = new TH1**[fNcentBins];
131  fHistJetTracksPt = new TH1**[fNcentBins];
132  fhTrackingEfficiency = new TProfile*[fNcentBins];
133  fHistPtaVsJt = new TProfile**[fNcentBins];
134  fHistBgPtaVsJt = new TProfile**[fNcentBins];
135  //CentBinBorders = new Double_t[10];
136 
137 
138  for (Int_t i = 0; i < fNcentBins; i++) {
139  fHistJTPta[i] = 0;
140  fHistLogJTPta[i] = 0;
141  fHistJTPta_all[i] = 0;
142  fHistJTBg[i] = 0;
143  fHistLogJTBg[i] = 0;
144  fHistJTPtaNonInv[i] = 0;
145  fHistLogJTPtaNonInv[i] = 0;
146  fHistJTPta_allNonInv[i] = 0;
147  fHistJTBgNonInv[i] = 0;
148  fHistLogJTBgNonInv[i] = 0;
149  fHistBackgroundDone[i] = 0;
150  fHistTracksPt[i] = 0;
151  fHistTracksJt[i] = 0;
152  fHistTracksEta[i] = 0;
153  fHistClustersPt[i] = 0;
154  fHistLeadingJetPt[i] = 0;
155  fHistJetsPt[i] = 0;
156  fHistJetsCorrPt[i] = 0;
157  fHistJetsCorrPtVsNonCorr[i] = 0;
158  fHistBgMulti[i] = 0;
159  fHistBgPt[i] = 0;
160  fHistJetEta[i] = 0;
161  fHistJetMulti[i] = 0;
162  fHistJetTracksPt[i] = 0;
163  fhTrackingEfficiency[i] = 0;
164  fHistPtaVsJt[i] = 0;
165  fHistBgPtaVsJt[i] = 0;
166  }
167 
168  /*for(Int_t i = 0; i < nJetsConts; i++){
169  fJetsConts[i] = 0;
170  }*/
171  SetMakeGeneralHistograms(kTRUE);
172 }
173 
174 //________________________________________________________________________
176  AliAnalysisTaskEmcalJet(name, kTRUE),
177  fHistTracksPt(0),
178  fHistTracksJt(0),
179  fHistTracksEta(0),
180  fHistClustersPt(0),
181  fHistLeadingJetPt(0),
182  fHistJetsPt(0),
183  fHistJetsCorrPt(0),
184  fHistJetsCorrPtVsNonCorr(0),
185  fHistBackgroundDone(0),
186  fHistJTPta(0),
187  fHistLogJTPta(0),
188  fHistJTPta_all(0),
189  fHistJTBg(0),
190  fHistLogJTBg(0),
191  fHistPtaVsJt(0),
192  fHistBgPtaVsJt(0),
193  fHistJTPtaNonInv(0),
194  fHistLogJTPtaNonInv(0),
195  fHistJTPta_allNonInv(0),
196  fHistJTBgNonInv(0),
197  fHistLogJTBgNonInv(0),
198  fHistBgMulti(0),
199  fHistBgPt(0),
200  fHistJetEta(0),
201  fHistJetMulti(0),
202  fHistJetTracksPt(0),
203  fhTrackingEfficiency(0),
204  fNpttBins(1),
205  fNptaBins(1),
206  fEffMode(1),
207  fJetsCont(0),
208  //fJetsConts(0),
209  //nJetsConts(0),
210  fTracksCont(0),
211  fCaloClustersCont(0),
212  fhVertexZ(0),
213  fHistEvtSelection(0),
214  fPrimaryVertex(0),
215  fTracks(0),
216  fTrackArrayName("nonejk"),
217  runPeriod(""),
218  fEfficiency(0),
219  fVertexHelper(0),
220  debug(0)
221 {
222  // Standard constructor.
223  fHistTracksPt = new TH1*[fNcentBins];
224  fHistTracksJt = new TH1*[fNcentBins];
225  fHistTracksEta = new TH1*[fNcentBins];
228  fHistJetsPt = new TH1**[fNcentBins];
229  fHistJetsCorrPt = new TH1**[fNcentBins];
230  fHistJetsCorrPtVsNonCorr = new TProfile*[fNcentBins];
232  fHistJTPta = new TH1***[fNcentBins];
233  fHistLogJTPta = new TH1***[fNcentBins];
234  fHistJTPta_all = new TH1***[fNcentBins];
235  fHistJTBg = new TH1***[fNcentBins];
236  fHistLogJTBg = new TH1***[fNcentBins];
237  fHistJTPtaNonInv = new TH1***[fNcentBins];
240  fHistJTBgNonInv = new TH1***[fNcentBins];
242  fHistBgMulti = new TH1**[fNcentBins];
243  fHistBgPt = new TH1**[fNcentBins];
244  fHistJetEta = new TH1**[fNcentBins];
245  fHistJetMulti = new TH1**[fNcentBins];
247  fhTrackingEfficiency = new TProfile*[fNcentBins];
248  fHistPtaVsJt = new TProfile**[fNcentBins];
249  fHistBgPtaVsJt = new TProfile**[fNcentBins];
250 
251 
252  for (Int_t i = 0; i < fNcentBins; i++) {
253  fHistJTPta[i] = 0;
254  fHistLogJTPta[i] = 0;
255  fHistJTPta_all[i] = 0;
256  fHistJTBg[i] = 0;
257  fHistLogJTBg[i] = 0;
258  fHistJTPtaNonInv[i] = 0;
259  fHistLogJTPtaNonInv[i] = 0;
260  fHistJTPta_allNonInv[i] = 0;
261  fHistJTBgNonInv[i] = 0;
262  fHistLogJTBgNonInv[i] = 0;
263  fHistBackgroundDone[i] = 0;
264  fHistTracksPt[i] = 0;
265  fHistTracksJt[i] = 0;
266  fHistTracksEta[i] = 0;
267  fHistClustersPt[i] = 0;
268  fHistLeadingJetPt[i] = 0;
269  fHistJetsPt[i] = 0;
270  fHistJetsCorrPt[i] = 0;
272  fHistBgMulti[i] = 0;
273  fHistBgPt[i] = 0;
274  fHistJetEta[i] = 0;
275  fHistJetMulti[i] = 0;
276  fHistJetTracksPt[i] = 0;
277  fhTrackingEfficiency[i] = 0;
278  fHistPtaVsJt[i] = 0;
279  fHistBgPtaVsJt[i] = 0;
280  }
281 
283 }
284 
285 //________________________________________________________________________
287 {
288  // Destructor.
289 }
290 
291 
293  fNcentBins=n;
294  if(debug > 0){
295  cout << "AliAnalysisTaskJetJTJT::setCentBinBorders: " << endl;
296  }
297  for(int i= 0 ; i < fNcentBins; i++){
298  CentBinBorders[i]= c[i];
299  if(debug > 0)
300  cout << CentBinBorders[i] << endl;
301  }
302 }
303 
305  fNpttBins=n;
306  if(debug > 0)
307  cout << "AliAnalysisTaskJetJTJT::setTriggPtBorders: " << endl;
308  for(int i= 0 ; i < fNpttBins; i++){
309  TriggPtBorders[i]= c[i];
310  if(debug > 0)
311  cout << TriggPtBorders[i] << endl;
312  }
313 }
314 
316  fNptaBins=n;
317  if(debug > 0)
318  cout << "AliAnalysisTaskJetJTJT::setAssocPtBorders: " << endl;
319  for(int i= 0 ; i < fNptaBins; i++){
320  AssocPtBorders[i]= c[i];
321  if(debug > 0)
322  cout << AssocPtBorders[i] << endl;
323  }
324 }
325 
326 
327 //________________________________________________________________________
329 {
330  // Create user output.
331 
333  if(debug > 0)
334  cout << "Creating Histograms" << endl;
335 
337  /*for(int i = 0; i < nJetsConts ; i++){
338  fJetsConts[0] = GetJetContainer(i);
339  }*/
340  /*if(fJetsCont) { //get particles and clusters connected to jets
341  fTracksCont = fJetsCont->GetParticleContainer();
342  fCaloClustersCont = fJetsCont->GetClusterContainer();
343  } else { //no jets, just analysis tracks and clusters
344  fTracksCont = GetParticleContainer(0);
345  fCaloClustersCont = GetClusterContainer(0);
346  }*/
349  fTracksCont->SetClassName("AliVTrack");
350  fCaloClustersCont->SetClassName("AliAODCaloCluster");
351 
352  TString histname;
353  //Int_t fMinBinJt = 0;
354  //Int_t fMaxBinJt = 5;
355 
356  Int_t NBINSJt=160;
357  double LogBinsJt[NBINSJt+1], LimLJt=0.001, LimHJt=10;
358  double logBWJt = (TMath::Log(LimHJt)-TMath::Log(LimLJt))/(NBINSJt-1);
359  LogBinsJt[0] = 0;
360  for(int ij=1;ij<=NBINSJt;ij++) LogBinsJt[ij]=LimLJt*exp(ij*logBWJt);
361 
362  Int_t NBINSPt = 160;
363  double LogBinsPt[NBINSPt+1], LimLPt=0.1, LimHPt=20;
364  double logBWPt = (TMath::Log(LimHPt)-TMath::Log(LimLPt))/(NBINSPt-1);
365  LogBinsPt[0] = 0;
366  for(int ib=1;ib<=NBINSPt;ib++) LogBinsPt[ib]=LimLPt*exp(ib*logBWJt);
367 
368  int NBINSJtW=160;
369  double LimLJtW=TMath::Log(LimLJt), LimHJtW=TMath::Log(LimHJt);
370 
371  //==== Efficiency ====
372  if(debug > 0)
373  cout << "AliAnalysisTaskJetJTJT::UserCreateOutputObjects: Creating efficiency" << endl;
375  fEfficiency->SetMode( fEffMode ); // 0:NoEff, 1:Period 2:RunNum 3:Auto
376  fEfficiency->SetDataPath("alien:///alice/cern.ch/user/d/djkim/legotrain/efficieny/data"); // Efficiency root file location local or alien
377 
378 
379  histname = "fHistVertexZ";
380  fhVertexZ = new TH1F(histname.Data(), histname.Data(), 200, -20., 20.);
381  fhVertexZ->GetXaxis()->SetTitle("#Delta z(cm)");
382  fhVertexZ->GetYaxis()->SetTitle("dN^{Events}/dz");
383  fOutput->Add(fhVertexZ);
384 
385  // Event statistics
386  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 4, -0.5, 4.5);
387  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
388  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
389  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"pile up (rejected)");
390  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
391  //fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
392 
394  for (Int_t ic = 0; ic < fNcentBins; ic++) {
395  if (fParticleCollArray.GetEntriesFast()>0) {
396  histname = "fHistTracksPt_";
397  histname += ic;
398  fHistTracksPt[ic] = new TH1F(histname.Data(), histname.Data(), NBINSPt, LogBinsPt);
399  fHistTracksPt[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
400  fHistTracksPt[ic]->GetYaxis()->SetTitle("tracks");
401  fOutput->Add(fHistTracksPt[ic]);
402 
403  histname = "fHistTracksJt_";
404  histname += ic;
405  fHistTracksJt[ic] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
406  fHistTracksJt[ic]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
407  fHistTracksJt[ic]->GetYaxis()->SetTitle("tracks");
408  fOutput->Add(fHistTracksJt[ic]);
409 
410  histname = "fHistTracksEta_";
411  histname += ic;
412  fHistTracksEta[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, -2, 2);
413  fHistTracksEta[ic]->GetXaxis()->SetTitle("#eta");
414  fHistTracksEta[ic]->GetYaxis()->SetTitle("tracks");
415  fOutput->Add(fHistTracksEta[ic]);
416  }
417 
418  histname = "fhTrackingEfficiency_";
419  histname += ic;
420  fhTrackingEfficiency[ic] = new TProfile(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
421  fhTrackingEfficiency[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
422  fhTrackingEfficiency[ic]->GetYaxis()->SetTitle("counts");
423  fOutput->Add(fhTrackingEfficiency[ic]);
424  fHistJTPta[ic] = new TH1**[fNpttBins];
425  fHistLogJTPta[ic] = new TH1**[fNpttBins];
426  fHistJTPta_all[ic] = new TH1**[fNpttBins];
427  fHistJTBg[ic] = new TH1**[fNpttBins];
428  fHistLogJTBg[ic] = new TH1**[fNpttBins];
429  fHistJTPtaNonInv[ic] = new TH1**[fNpttBins];
430  fHistLogJTPtaNonInv[ic] = new TH1**[fNpttBins];
431  fHistJTPta_allNonInv[ic] = new TH1**[fNpttBins];
432  fHistJTBgNonInv[ic] = new TH1**[fNpttBins];
433  fHistLogJTBgNonInv[ic] = new TH1**[fNpttBins];
434  fHistJetsPt[ic] = new TH1*[fNpttBins];
435  fHistJetsCorrPt[ic] = new TH1*[fNpttBins];
436  fHistBackgroundDone[ic] = new TH1*[fNpttBins];
437  fHistBgMulti[ic] = new TH1*[fNpttBins];
438  fHistBgPt[ic] = new TH1*[fNpttBins];
439  fHistJetEta[ic] = new TH1*[fNpttBins];
440  fHistJetMulti[ic] = new TH1*[fNpttBins];
441  fHistJetTracksPt[ic] = new TH1*[fNpttBins];
442  fHistPtaVsJt[ic] = new TProfile*[fNpttBins];
443  fHistBgPtaVsJt[ic] = new TProfile*[fNpttBins];
444  for(Int_t j=0; j < fNpttBins; j++){
445  fHistJTPta[ic][j] = new TH1*[fNptaBins];
446  fHistLogJTPta[ic][j] = new TH1*[fNptaBins];
447  fHistJTPta_all[ic][j] = new TH1*[fNptaBins];
448  fHistJTBg[ic][j] = new TH1*[fNptaBins];
449  fHistLogJTBg[ic][j] = new TH1*[fNptaBins];
450  fHistJTPtaNonInv[ic][j] = new TH1*[fNptaBins];
451  fHistLogJTPtaNonInv[ic][j] = new TH1*[fNptaBins];
452  fHistJTPta_allNonInv[ic][j] = new TH1*[fNptaBins];
453  fHistJTBgNonInv[ic][j] = new TH1*[fNptaBins];
454  fHistLogJTBgNonInv[ic][j] = new TH1*[fNptaBins];
455  for(Int_t k=0; k < fNptaBins; k++){
456  fHistJTPta[ic][j][k] = 0;
457  fHistLogJTPta[ic][j][k] = 0;
458  fHistJTPta_all[ic][j][k] = 0;
459  fHistJTBg[ic][j][k] = 0;
460  fHistLogJTBg[ic][j][k] = 0;
461  fHistJTPtaNonInv[ic][j][k] = 0;
462  fHistLogJTPtaNonInv[ic][j][k] = 0;
463  fHistJTPta_allNonInv[ic][j][k] = 0;
464  fHistJTBgNonInv[ic][j][k] = 0;
465  fHistLogJTBgNonInv[ic][j][k] = 0;
466  }
467  fHistJetsPt[ic][j] = 0;
468  fHistJetsCorrPt[ic][j] = 0;
469  fHistBackgroundDone[ic][j] = 0;
470  fHistBgMulti[ic][j] = 0;
471  fHistBgPt[ic][j] = 0;
472  fHistJetEta[ic][j] = 0;
473  fHistJetMulti[ic][j] =0;
474  fHistJetTracksPt[ic][j] = 0;
475  fHistPtaVsJt[ic][j] = 0;
476  fHistBgPtaVsJt[ic][j] = 0;
477  }
478 
479 
480  if (fParticleCollArray.GetEntriesFast()>0) {
481  for(Int_t iptt = 0 ; iptt < fNpttBins; iptt++){
482  for(Int_t ipta = 0 ; ipta < fNptaBins; ipta++){
483  histname = "hJTPtaD00C";
484  //histname += ic;
485  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
486  if(debug > 1)
487  cout << histname << endl;
488  fHistJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
489  fHistJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
490  fHistJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
491  fOutput->Add(fHistJTPta[ic][iptt][ipta]);
492 
493  histname = "hLogJTPtaD00C";
494  //histname += ic;
495  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
496  if(debug > 1)
497  cout << histname << endl;
498  fHistLogJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
499  fHistLogJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
500  fHistLogJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
501  fOutput->Add(fHistLogJTPta[ic][iptt][ipta]);
502 
503  histname = "hJTPta_allD00C";
504  //histname += ic;
505  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
506  if(debug > 1)
507  cout << histname << endl;
508  fHistJTPta_all[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
509  fHistJTPta_all[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
510  fHistJTPta_all[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
511  fOutput->Add(fHistJTPta_all[ic][iptt][ipta]);
512 
513  histname = "hJTBg";
514  //histname += ic;
515  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
516  if(debug > 1)
517  cout << histname << endl;
518  fHistJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
519  fHistJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
520  fHistJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
521  fOutput->Add(fHistJTBg[ic][iptt][ipta]);
522 
523  histname = "hLogJTBg";
524  //histname += ic;
525  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
526  if(debug > 1)
527  cout << histname << endl;
528  fHistLogJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
529  fHistLogJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
530  fHistLogJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
531  fOutput->Add(fHistLogJTBg[ic][iptt][ipta]);
532 
533  histname = "hJTPtaNonInvD00C";
534  //histname += ic;
535  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
536  if(debug > 1)
537  cout << histname << endl;
538  fHistJTPtaNonInv[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
539  fHistJTPtaNonInv[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
540  fHistJTPtaNonInv[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
541  fOutput->Add(fHistJTPtaNonInv[ic][iptt][ipta]);
542 
543  histname = "hLogJTPtaNonInvD00C";
544  //histname += ic;
545  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
546  if(debug > 1)
547  cout << histname << endl;
548  fHistLogJTPtaNonInv[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
549  fHistLogJTPtaNonInv[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
550  fHistLogJTPtaNonInv[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
551  fOutput->Add(fHistLogJTPtaNonInv[ic][iptt][ipta]);
552 
553  histname = "hJTPta_allNonInvD00C";
554  //histname += ic;
555  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
556  if(debug > 1)
557  cout << histname << endl;
558  fHistJTPta_allNonInv[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
559  fHistJTPta_allNonInv[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
560  fHistJTPta_allNonInv[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
561  fOutput->Add(fHistJTPta_allNonInv[ic][iptt][ipta]);
562 
563  histname = "hJTBgNonInv";
564  //histname += ic;
565  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
566  if(debug > 1)
567  cout << histname << endl;
568  fHistJTBgNonInv[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
569  fHistJTBgNonInv[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
570  fHistJTBgNonInv[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
571  fOutput->Add(fHistJTBgNonInv[ic][iptt][ipta]);
572 
573  histname = "hLogJTBg";
574  //histname += ic;
575  histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
576  if(debug > 1)
577  cout << histname << endl;
578  fHistLogJTBgNonInv[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
579  fHistLogJTBgNonInv[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
580  fHistLogJTBgNonInv[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
581  fOutput->Add(fHistLogJTBgNonInv[ic][iptt][ipta]);
582 
583 
584  }
585  }
586  }
587 
588  if (fClusterCollArray.GetEntriesFast()>0) {
589  histname = "fHistClustersPt_";
590  histname += ic;
591  fHistClustersPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
592  fHistClustersPt[ic]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
593  fHistClustersPt[ic]->GetYaxis()->SetTitle("counts");
594  fOutput->Add(fHistClustersPt[ic]);
595  }
596 
597  if (fJetCollArray.GetEntriesFast()>0) {
598  histname = "fHistLeadingJetPt_";
599  histname += ic;
600  fHistLeadingJetPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
601  fHistLeadingJetPt[ic]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
602  fHistLeadingJetPt[ic]->GetYaxis()->SetTitle("counts");
603  fOutput->Add(fHistLeadingJetPt[ic]);
604 
605  histname = "fHistJetsCorrPtVsNonCorr_";
606  histname += Form("C%02d", ic);
607  fHistJetsCorrPtVsNonCorr[ic] = new TProfile(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
608  fHistJetsCorrPtVsNonCorr[ic]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
609  fHistJetsCorrPtVsNonCorr[ic]->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
611 
612  for(Int_t iptt = 0 ; iptt < fNpttBins; iptt++){
613 
614  histname = "fHistJetsPt_";
615  histname += Form("C%02dT%02d", ic, iptt);
616  fHistJetsPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
617  fHistJetsPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
618  fHistJetsPt[ic][iptt]->GetYaxis()->SetTitle("counts");
619  fOutput->Add(fHistJetsPt[ic][iptt]);
620 
621  histname = "fHistJetsCorrPt_";
622  histname += Form("C%02dT%02d", ic, iptt);
623  fHistJetsCorrPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
624  fHistJetsCorrPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
625  fHistJetsCorrPt[ic][iptt]->GetYaxis()->SetTitle("counts");
626  fOutput->Add(fHistJetsCorrPt[ic][iptt]);
627 
628  histname = "fHistBackgroundDone_";
629  histname += Form("C%02dT%02d", ic, iptt);;
630  fHistBackgroundDone[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 2, -1, 2);
631  fHistBackgroundDone[ic][iptt]->GetXaxis()->SetTitle("Number of jets");
632  fHistBackgroundDone[ic][iptt]->GetYaxis()->SetTitle("0 = not done, 1 = done");
633  fOutput->Add(fHistBackgroundDone[ic][iptt]);
634 
635  histname = "fHistJetEta_";
636  histname += Form("C%02dT%02d", ic, iptt);
637  fHistJetEta[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, -2, 2);
638  fHistJetEta[ic][iptt]->GetXaxis()->SetTitle("#eta");
639  fHistJetEta[ic][iptt]->GetYaxis()->SetTitle("jets");
640  fOutput->Add(fHistJetEta[ic][iptt]);
641 
642  histname = "fHistJetMulti_";
643  histname += Form("C%02dT%02d", ic, iptt);
644  fHistJetMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
645  fHistJetMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
646  fHistJetMulti[ic][iptt]->GetYaxis()->SetTitle("jets");
647  fOutput->Add(fHistJetMulti[ic][iptt]);
648 
649  histname = "fHistBgMulti_";
650  histname += Form("C%02dT%02d", ic, iptt);
651  fHistBgMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
652  fHistBgMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
653  fHistBgMulti[ic][iptt]->GetYaxis()->SetTitle("Events");
654  fOutput->Add(fHistBgMulti[ic][iptt]);
655 
656  histname = "fHistBgPt_";
657  histname += Form("C%02dT%02d", ic, iptt);
658  fHistBgPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins,fMinBinPt, fMaxBinPt/4);
659  fHistBgPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}");
660  fHistBgPt[ic][iptt]->GetYaxis()->SetTitle("tracks");
661  fOutput->Add(fHistBgPt[ic][iptt]);
662 
663  histname = "fHistJetTracksPt_";
664  //histname += ic;
665  histname += Form("C%02dT%02d", ic, iptt);
666  if(debug > 1)
667  cout << histname << endl;
668  fHistJetTracksPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt/10);
669  fHistJetTracksPt[ic][iptt]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
670  fHistJetTracksPt[ic][iptt]->GetYaxis()->SetTitle("counts");
671  fOutput->Add(fHistJetTracksPt[ic][iptt]);
672 
673  histname = "fHistPtaVsJt_";
674  histname += Form("C%02dT%02d", ic, iptt);
675  fHistPtaVsJt[ic][iptt] = new TProfile(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
676  fHistPtaVsJt[ic][iptt]->GetXaxis()->SetTitle("j_{T,track} (GeV/c)");
677  fHistPtaVsJt[ic][iptt]->GetYaxis()->SetTitle("p_{T,track} (GeV/c)");
678  fOutput->Add(fHistPtaVsJt[ic][iptt]);
679 
680  histname = "fHistBgPtaVsJt_";
681  histname += Form("C%02dT%02d", ic, iptt);
682  fHistBgPtaVsJt[ic][iptt] = new TProfile(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
683  fHistBgPtaVsJt[ic][iptt]->GetXaxis()->SetTitle("j_{T,track} (GeV/c)");
684  fHistBgPtaVsJt[ic][iptt]->GetYaxis()->SetTitle("p_{T,track} (GeV/c)");
685  fOutput->Add(fHistBgPtaVsJt[ic][iptt]);
686  }
687 
688  /*
689  if (!(GetJetContainer()->GetRhoName().IsNull())) {
690  histname = "fHistJetsCorrPtArea_";
691  histname += i;
692  fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
693  fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
694  fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
695  fOutput->Add(fHistJetsCorrPtArea[i]);
696  }
697  */
698  }
699  }
700 
701  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
702 }
703 
704 //________________________________________________________________________
706 {
707  // Fill histograms.
708 
709  fHistEvtSelection->Fill(0); //Count accepted events
710  AliCentrality *aliCent = InputEvent()->GetCentrality();
711  fCentBin = 0;
712  if (aliCent) {
713  //fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
714  fCent = aliCent->GetCentralityPercentile("V0M");
715  /*if(debug > 0){
716  cout << "Centrality " << fCent << endl;
717  }*/
718  for(int ic = 0; ic < fNcentBins; ic++){
719  /*if(debug > 0){
720  cout << "ic: " << ic << " / " << fNcentBins << endl;
721  cout << "Centrality bin " << fCentBin << endl;
722  cout << "Border: " << CentBinBorders[ic] << endl;
723  } */
724  if(fCent > CentBinBorders[ic]){
725  fCentBin = ic;
726  }
727  }
728  //cout << "Centrality bin: " << fCentBin << endl;
729  } else {
730  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
731  fCentBin = 3;
732  }
733  int fHadronSelectionCut = 5; //5=Hybrid cut
734  if (fTracksCont) {
735  fTracksCont->ResetCurrentID();
736  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
737  while(track) {
738  double ptt = track->Pt();
739 
740  //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
741  //double effCorr = 1./fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fCent); // here you generate warning if ptt>30
742  if(debug > 0)
743  cout << "Getting efficiency correction for ptt " << ptt << " with centrality " << fCent << endl;
744  double effCorr = fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fCent); // here you generate warning if ptt>30
745  //double effCorr = 1.;
746  if(debug > 0)
747  cout << "Filling fhTrackingEfficiency fCentBin: " << fCentBin << " ptt: " << ptt << " with efficiency: " << effCorr << endl;
748  fhTrackingEfficiency[fCentBin]->Fill( ptt, effCorr );
749  //triggTr->SetTrackEff( 1./effCorr );
750  //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
751 
752  if(ptt > 0 && 1.0/ptt > 0){
753  fHistTracksPt[fCentBin]->Fill(ptt,effCorr);
754  fHistTracksEta[fCentBin]->Fill(track->Eta(),effCorr);
755  }
756 
757 
758  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
759  }
760  }
761 
762  if (fCaloClustersCont) {
763  fCaloClustersCont->ResetCurrentID();
764  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
765  while(cluster) {
766  TLorentzVector nPart;
767  cluster->GetMomentum(nPart, fVertex);
768  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
769 
771  }
772  }
773 
774  Int_t fPttBin, fPtaBin;
775  fPtaBin = 0;
776 
777 
778  if (fJetsCont) {
779  //Int_t Njets = fJetsCont->GetNJets();
780  Int_t Njets = 0;
781  if(debug > 1){
782  cout << "Number of Jets: " << Njets << endl;
783  }
784 
785  //Make arrays to hold jets to be tested in background jT
786  Float_t jetPhis[200] = {};
787  Float_t jetEtas[200] = {};
788  fJetsCont->ResetCurrentID();
790  Int_t ij = 0;
791  while(jet) {
792  //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
793  if(jet->Pt() > 5){ //Only consider jets with pT > 5 GeV
794  jetPhis[ij] = jet->Phi();
795  jetEtas[ij] = jet->Eta();
796  ij++;
797  Njets++;
798  if(debug > 1)
799  cout << "i: " << ij << " jetPhi: " << jetPhis[ij] << " jetEta: " << jetEtas[ij] << endl;
800  }else{
801  //jetPhis[ij] = 100;
802  //jetEtas[ij] = 100;
803  //Njets--;
804  if(debug > 1)
805  cout << "jetPt: " << jet->Pt() << " jetPhi: " << jet->Phi() << " jetEta: " << jet->Eta() << endl;
806  }
807  //i++;
808  jet = fJetsCont->GetNextAcceptJet();
809  }
810 
811 
812  //fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
813  fJetsCont->ResetCurrentID();
814  jet = fJetsCont->GetNextAcceptJet();
815  while(jet) {
816  if(jet->Pt() > 5){
817  if(jet->Eta() < -0.4 || jet->Eta() > 0.4){ //TODO Fix
818  if(debug > 0)
819  cout << "Jet outside eta range, Eta: " << jet->Eta() << endl;
820  jet = fJetsCont->GetNextAcceptJet();
821  continue;
822  }
823  //Get the trigger pT bin
824  fPttBin = 0;
825  for(int iptt = 0 ; iptt < fNpttBins; iptt++){
826  if(jet->Pt() > TriggPtBorders[iptt]){
827  fPttBin = iptt;
828  }
829 
830  }
831  fHistJetEta[fCentBin][fPttBin]->Fill(jet->Eta());
832  Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
833  if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
834  fHistJetsPt[fCentBin][fPttBin]->Fill(jet->Pt(),1.0/jet->Pt()); //Fill jet dN/(pT dpT)
835  fHistJetsCorrPt[fCentBin][fPttBin]->Fill(corrPt,1.0/corrPt); //Fill jet dN/(pT dpT)
836  fHistJetsCorrPtVsNonCorr[fCentBin]->Fill(jet->Pt(),corrPt); //Fill jet dN/(pT dpT)
837 
838  Int_t nTrack = jet->GetNumberOfTracks();
839  if (debug > 0)
840  cout << "Number of tracks " << nTrack << " Jet Pt: " << jet->Pt() << endl;
841  fHistJetMulti[fCentBin][fPttBin]->Fill(nTrack);
842  for(Int_t it = 0; it < nTrack; it++ ){
843  AliVParticle *track = (AliVParticle*)jet->TrackAt( it, fTracks );
844  if( !track ){
845  cout << "No Track found" << endl;
846  continue;
847  }
848  fPtaBin = 0; //Get the associated pT bin
849  for(int ipta = 0 ; ipta < fNptaBins; ipta++){
850  if(track->Pt() > AssocPtBorders[ipta]){
851  fPtaBin = ipta;
852  }
853  }
854  fHistJetTracksPt[fCentBin][fPttBin]->Fill(track->Pt());
855  if(debug > 2)
856  cout << "Filling fHistJetTracksPt C" << fCentBin << " T" << fPttBin << endl;
857  Float_t jt = getJt(track,jet,0);
858  double effCorr = fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent); // here you generate warning if ptt>30
859  if(jt > 0 && 1.0/jt > 0){
860  fHistTracksJt[fCentBin]->Fill(jt,effCorr/jt); //Fill dN/(djT jT)
861  fHistJTPta[fCentBin][fPttBin][fPtaBin]->Fill(jt,effCorr/jt); //Fill dN/(djT jT)
862  fHistLogJTPta[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),effCorr/(jt*jt)); //Fill logarithmic dN/(dln(jT) jT^2)
863  fHistPtaVsJt[fCentBin][fPttBin]->Fill(jt,track->Pt(),effCorr); //Fill j_T vs p_Ta histogram
864  fHistJTPtaNonInv[fCentBin][fPttBin][fPtaBin]->Fill(jt,effCorr); //Fill dN/(djT)
865  fHistLogJTPtaNonInv[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0*effCorr); //Fill logarithmic dN/(dln(jT))
866  }
867  if(debug > 1)
868  cout << "Filling JT C" << fCentBin << "T" << fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << effCorr/jt<< endl;
869  }
870 
871  //Get Jet azimuth and rapidity of jet
872  Float_t jetAngle = jet->Phi();
873  Float_t jetRap = jet->Eta();
874 
875  //Rotate jet angle for background cone
876  Float_t rotatedAngle = jetAngle+TMath::Pi()/2;
877  if(rotatedAngle > TMath::Pi()*2){
878  rotatedAngle = rotatedAngle- TMath::Pi()*2;
879  }
880  AliEmcalJet *bgCone = new AliEmcalJet(jet->Pt(), jetRap, rotatedAngle, jet->M());
881 
882  Float_t jetArea = jet->Area();
883  Float_t testRadius = TMath::Sqrt(jetArea/TMath::Pi());
884 
885  Bool_t doBg = 1;
886 
887  //Test if there are jets in the background test cone
888  for(int i_j = 0; i_j < Njets; i_j++){
889  //Float_t diffR = TMath::Sqrt(TMath::Power(jetPhis[i_j]-rotatedAngle,2)+TMath::Power(jetEtas[i_j]-jetRap,2));
890  Float_t diffR = getDiffR(jetPhis[i_j],rotatedAngle,jetEtas[i_j],jetRap);
891  if(debug > 1){
892  cout << "i_j: " << i_j << " JetPhi: " << jetPhis[i_j] << " jetEta: " << jetEtas[i_j] << endl;
893  cout << "DiffR: " << diffR << " doBG: " << doBg <<endl;
894  }
895  if(diffR < testRadius *2){ //Jets muts be at least 2*cone radius away from the background cone axis
896  doBg =0;
897  break;
898  }
899 
900  }
901 
902  // Do jT for all particles in respect to jet axis
903  if (fTracksCont) {
904  int counter = 0;
905  fTracksCont->ResetCurrentID();
906  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
907  while(track) {
908  Double_t jt = getJt(track,bgCone,0);
909  double effCorr = fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent); // here you generate warning if ptt>30
910  if(jt > 0 && 1.0/jt > 0){
911  fHistJTPta_all[fCentBin][fPttBin][fPtaBin]->Fill(jt,effCorr/jt);
912  fHistJTPta_allNonInv[fCentBin][fPttBin][fPtaBin]->Fill(jt,effCorr);
913  }
914  for(int ipta = 0 ; ipta < fNptaBins; ipta++){
915  if(track->Pt() > AssocPtBorders[ipta]){
916  fPtaBin = ipta;
917  }
918  }
919  //If background is to be filled
920  if(doBg){
921  //Float_t diffR = TMath::Sqrt(TMath::Power(track->Phi()-rotatedAngle,2)+TMath::Power(track->Eta()-jetRap,2));
922  Float_t diffR = getDiffR(track->Phi(),rotatedAngle,track->Eta(),jetRap);
923  //Particles in the rotated cone
924  if(diffR < testRadius){
925  counter++;
926  fHistBgPt[fCentBin][fPttBin]->Fill(track->Pt(),effCorr);
927  jt = getJt(track,bgCone,0);
928  if(jt > 0 && 1.0/jt > 0){
929  fHistJTBg[fCentBin][fPttBin][fPtaBin]->Fill(jt,effCorr/jt);
930  fHistLogJTBg[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),effCorr/(jt*jt));
931  fHistJTBgNonInv[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0*effCorr);
932  fHistLogJTBgNonInv[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0*effCorr);
933  fHistBgPtaVsJt[fCentBin][fPttBin]->Fill(jt,track->Pt(),effCorr);
934  }
935  if(debug > 1)
936  cout << "Filling Background C" << fCentBin << "T" << fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << effCorr/jt<< endl;
937  //Fill background jT
938  }
939  }
940  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
941  }
942  if(doBg){
943  fHistBgMulti[fCentBin][fPttBin]->Fill(counter);
944  }
945  }
946  if(doBg){
947  fHistBackgroundDone[fCentBin][fPttBin]->Fill(1);
948  }else{
949  fHistBackgroundDone[fCentBin][fPttBin]->Fill(0);
950  }
951  }
952  }
953  jet = fJetsCont->GetNextAcceptJet();
954  }
955  jet = fJetsCont->GetLeadingJet();
956  if(jet){
957  if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
958  fHistLeadingJetPt[fCentBin]->Fill(jet->Pt(),1.0/jet->Pt());
959  }
960  }
961  }
963  return kTRUE;
964 }
965 
966 
967 //-----------------------------------------------------------------------
968 Double_t AliAnalysisTaskJetJTJT::getJt(AliVTrack *track, AliEmcalJet *jet,int reverse){
969  Float_t dotproduct = 0;
970  Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
971  if(reverse){
972  dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
973  } else{
974  dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
975  }
976  Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
977  Float_t normproduct = constp*jetp;
978  Float_t costheta2 = dotproduct/normproduct;
979  //Float_t sintheta = sqrt(1-costheta2*costheta2);
980  Float_t jt = constp*sqrt(1-costheta2*costheta2);
981  return jt;
982 }
983 
984 Double_t AliAnalysisTaskJetJTJT::getJt(AliVParticle *track, AliEmcalJet *jet,int reverse){
985  Float_t dotproduct = 0;
986  Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
987  if(reverse){
988  dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
989  } else{
990  dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
991  }
992  Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
993  Float_t normproduct = constp*jetp;
994  Float_t costheta2 = dotproduct/normproduct;
995  //Float_t sintheta = sqrt(1-costheta2*costheta2);
996  Float_t jt = constp*sqrt(1-costheta2*costheta2);
997  return jt;
998 }
999 
1000 //Phi1 and Phi2 between 0 and 2 pi
1001 Double_t AliAnalysisTaskJetJTJT::getDiffR(double phi1, double phi2, double eta1, double eta2){
1002  Double_t diffPhi = TMath::Abs(phi1-phi2);
1003  if(diffPhi > TMath::Pi()){
1004  diffPhi = 2*TMath::Pi() - diffPhi;
1005  }
1006  return TMath::Sqrt(TMath::Power(diffPhi,2)+TMath::Power(eta1-eta2,2));
1007 }
1008 
1009 //________________________________________________________________________
1011 {
1012 
1014  return;
1015 
1016  Double_t deta = 999;
1017  Double_t dphi = 999;
1018 
1019  //Get closest cluster to track
1020  fTracksCont->ResetCurrentID();
1021  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
1022  while(track) {
1023  //Get matched cluster
1024  Int_t emc1 = track->GetEMCALcluster();
1025  if(fCaloClustersCont && emc1>=0) {
1026  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
1027  if(clusMatch) {
1028  AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
1029  //fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
1030  }
1031  }
1032  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
1033  }
1034 
1035  //Get closest track to cluster
1036  fCaloClustersCont->ResetCurrentID();
1037  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
1038  while(cluster) {
1039  TLorentzVector nPart;
1040  cluster->GetMomentum(nPart, fVertex);
1041  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
1042 
1043  //Get matched track
1044  AliVTrack *mt = NULL;
1045  AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
1046  if(acl) {
1047  if(acl->GetNTracksMatched()>1)
1048  mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
1049  }
1050  else {
1051  AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
1052  Int_t im = ecl->GetTrackMatchedIndex();
1053  if(fTracksCont && im>=0) {
1054  mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
1055  }
1056  }
1057  if(mt) {
1058  AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
1059  //fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
1060 
1061  //debugging
1062  /*
1063  if(mt->IsEMCAL()) {
1064  Int_t emc1 = mt->GetEMCALcluster();
1065  Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
1066  AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
1067  AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
1068  Printf("deta: %f dphi: %f",deta,dphi);
1069  }
1070  */
1071  }
1073  }
1074 }
1075 
1076 //________________________________________________________________________
1077 /*AliJetContainer* AliAnalysisTaskJetJTJT::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) {
1078 
1079  AliAnalysisTaskEmcalJet::ExecOnce();
1080  nJetsConts++;
1081  AliJetContainer *cont = 0x0;
1082  cont = AliAnalysisTaskEmcalJet::AddJetContainer(n,defaultCutType,jetRadius);
1083  return cont;
1084  }*/
1085 
1086 
1087 //________________________________________________________________________
1089 
1090  if(debug > 0){
1091  cout << "AliAnalysisTaskJetJTJT::ExecOnce(): " << endl;
1092  cout << "Get fTracks from " << fTrackArrayName.Data() << endl;
1093  }
1094  fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
1095 
1097  if(debug > 1)
1098  cout << "Efficiency: Set Run Period Name " << runPeriod << endl;
1100  if(debug > 1)
1101  cout << "Efficiency: Set Run number " << InputEvent()->GetRunNumber() << endl;
1102  fEfficiency->SetRunNumber( InputEvent()->GetRunNumber() ); //TODO Get run Number
1103  if(debug > 1)
1104  cout << "Efficiency: Load()" << endl;
1105  fEfficiency->Load();
1106  if(debug > 1)
1107  cout << "fEfficiency loaded" << endl;
1108 
1109 
1110 
1111  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
1112  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
1113  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
1114 
1115  // Initialize helper class (for vertex selection & pile up correction)
1116  fVertexHelper = new AliAnalysisUtils();
1117  fVertexHelper->SetCutOnZVertexSPD(kTRUE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
1118  fVertexHelper->SetMinVtxContr( 2 ); //Copied from Jiri
1119  fVertexHelper->SetMaxVtxZ( 10 ); //Copied from Jiri
1120 }
1121 
1122 //________________________________________________________________________
1124 {
1125  // Run analysis code here, if needed. It will be executed before FillHistograms().
1126  fHistEvtSelection->Fill(1); //Count input event
1127 
1128 
1129  //Vertex cut, z must be < 10cm
1130 
1131  if(!fVertexHelper || fVertexHelper->IsPileUpEvent(InputEvent())){
1132  fHistEvtSelection->Fill(2); //count events rejected by pileup
1133  return kFALSE;
1134  }
1135  fPrimaryVertex = InputEvent()->GetPrimaryVertex();
1136  if((TMath::Abs(fPrimaryVertex->GetZ()) > 10.0)){
1137  fHistEvtSelection->Fill(3); //count events rejected by vertex cut
1138  return kFALSE;
1139  }
1140 
1141  fhVertexZ->Fill(fPrimaryVertex->GetZ());
1142  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
1143 }
1144 
1145 //________________________________________________________________________
1147 {
1148  // Called once at the end of the analysis.
1149 }
1150 
1151 
1152 
1153 //________________________________________________________________________
1155  fMode(kAuto),
1156  fPeriod(-1),
1157  fDataPath(""),
1158  fName(""),
1159  fPeriodStr(""),
1160  fMCPeriodStr(""),
1161  fRunNumber(0),
1162  fTag(""),
1163  fInputRootName(""),
1164  fInputRoot(NULL),
1165  fCentBinAxis(0x0)
1166 {
1167  for (Int_t i=0; i < 3; i++) fEffDir[i] = 0;
1168 
1169 }
1170 
1172  fMode(obj.fMode),
1173  fPeriod(obj.fPeriod),
1174  fDataPath(obj.fDataPath),
1175  fName(obj.fName),
1176  fPeriodStr(obj.fPeriodStr),
1177  fMCPeriodStr(obj.fMCPeriodStr),
1178  fRunNumber(obj.fRunNumber),
1179  fTag(obj.fTag),
1180  fInputRootName(obj.fInputRootName),
1181  fInputRoot(obj.fInputRoot),
1182  fCentBinAxis(obj.fCentBinAxis)
1183 {
1184  // copy constructor TODO: handling of pointer members
1185  JUNUSED(obj);
1186  for (Int_t i=0; i < 3; i++) fEffDir[i] = obj.fEffDir[i];
1187 }
1188 
1190  // equal sign operator TODO: content
1191  JUNUSED(obj);
1192  return *this;
1193 }
1194 
1195 
1196 //________________________________________________________________________
1197 double JTJTEfficiency::GetCorrection( double pt, int icut , double cent ) const {
1198  if( fMode == kNotUse ) return 1;
1199  int icent = fCentBinAxis->FindBin( cent ) -1 ;
1200  if( icent < 0 || icent > fCentBinAxis->GetNbins()-1 ) {
1201  cout<<"J_WARNING : Centrality "<<cent<<" is out of CentBinBorder"<<endl;
1202  return 1;
1203  }
1204  // TODO error check for icent;
1205  int ivtx = 0;
1206  if( ! fCorrection[ivtx][icent][icut] ) {
1207  cout<<"J_WARNING : No Eff Info "<<pt<<"\t"<<icut<<"\t"<<cent<<"\t"<<icent<<endl;
1208  return 1;
1209  }
1210  TGraphErrors * gr = fCorrection[ivtx][icent][icut];
1211  //=== TEMPERORY SETTING. IT will be removed soon.
1212  if( pt > 30 ) pt = 30; // Getting eff of 30GeV for lager pt
1213  double cor = gr->Eval(pt);
1214  if ( cor < 0.2 ) cor = 0.2;
1215  return cor;
1216 }
1217 
1218 
1220  /*
1221  1. kNotUse : no Load, efficiency is 1 always
1222  2. has fInputRootName : Load that or crash
1223  3. has fName : Load fName [+runnumber] or crash
1224  4. has runnumber : Find Good MC period from AliJRunTable, or crash
1225  3. has period : Find Good MC period from AliJRunTable, or crash
1226 
1227  }
1228  */
1229 if(fPeriodStr == "LHC10b"){
1230  fInputRootName = "Eff--LHC10b-LHC10d1-0-.root";
1231 }
1232 if(fPeriodStr == "LHC10c"){
1233  fInputRootName = "Eff--LHC10c-LHC10d4-0-.root";
1234 }
1235 if(fPeriodStr == "LHC10d"){
1236  fInputRootName = "Eff--LHC10d-LHC10f6a-0-.root";
1237 }
1238 if(fPeriodStr == "LHC10e"){
1239  fInputRootName = "Eff--LHC10e-LHC10e20-0-.root";
1240 }
1241 if(fPeriodStr == "LHC10h"){
1242  fInputRootName = "Eff--LHC10h-LHC11a10a_bis-0-.root";
1243 }
1244 if(fPeriodStr == "LHC11a"){
1245  fInputRootName = "Eff--LHC11a-LHC11b10a-0-.root";
1246 }
1247 if(fPeriodStr == "LHC13b"){
1248  fInputRootName = "Eff--LHC13b-LHC13b2-efix_p1-0-.root";
1249 }
1250 
1251 if(fPeriodStr == "LHC13c"){
1252  fInputRootName = "Eff--LHC13c-LHC13b2-efix_p1-0-.root";
1253 }
1254 if(fPeriodStr == "LHC13d"){
1255  fInputRootName = "Eff--LHC13d-LHC13b2-efix_p1-0-.root";
1256 }
1257 if(fPeriodStr == "LHC13e"){
1258  fInputRootName = "Eff--LHC13e-LHC13b2-efix_p1-0-.root";
1259 }
1260 
1261 return fInputRootName;
1262 }
1263 
1265  GetEffName();
1267  return fInputRootName;
1268 }
1269 
1270 
1271 //________________________________________________________________________
1273  // Load Efficiency File based on fMode
1274  if( fMode == kNotUse ) {
1275  cout<<"J_WARNING : Eff Mode is \"NOTUSE\". eff is 1 !!!"<<endl;
1276  return true;
1277  }
1278  GetEffFullName();
1279  if (TString(fInputRootName).BeginsWith("alien:")) TGrid::Connect("alien:");
1280  fInputRoot = TFile::Open( fInputRootName);
1281  //fInputRoot = new TFile( fInputRootName,"READ");
1282  if( !fInputRoot ) {
1283  cout << "J_ERROR : %s does not exist" << fInputRootName << endl;
1284  return false;
1285  }
1286 
1287  //fEffDir[0] = (TDirectory*)fInputRoot->Get("EffRE");
1289  fEffDir[2] = (TDirectory*)fInputRoot->Get("Efficiency");
1290  //iif( fEffDir[0] && fEffDir[1] && fEffDir[2] )
1291  if( !fEffDir[2] )
1292  {
1293  cout << "J_ERROR : Directory EFF is not exist"<<endl;
1294  return false;
1295  }
1296 
1297  fCentBinAxis = (TAxis*)fEffDir[2]->Get("CentralityBin");
1298  if( !fCentBinAxis ){
1299  cout << "J_ERROR : No CentralityBin in directory" << endl;
1300  return false;
1301  }
1302 
1303 
1304  int nVtx = 1;
1305  int nCentBin = fCentBinAxis->GetNbins();
1306  for( int ivtx=0;ivtx<nVtx;ivtx++ ){
1307  for( int icent=0;icent<nCentBin;icent++ ){
1308  for( int icut=0;icut<kJNTrackCuts;icut++ ){
1309  fCorrection[ivtx][icent][icut]
1310  = (TGraphErrors*) fEffDir[2]->Get(Form("gCor%02d%02d%02d", ivtx,icent,icut));
1311  //cout<<"J_LOG : Eff graph - "<<Form("gCor%02d%02d%02d", ivtx,icent,icut)<<" - "<<g<<endl;
1312  }
1313  }
1314  }
1315  cout<<"J_LOG : Eff file is "<<fInputRootName<<endl;
1316  cout<<"J_LOG : Eff Cent Bins are ";
1317  for( int i=0;i<=nCentBin;i++ ){
1318  //cout<<fCentBinAxis->GetXbins()->At(i)<<" ";
1319  }
1320  //cout<<endl;
1321  return true;
1322 }
1323 
1324 
1325 
TH1 **** fHistJTPta_all
Logarithmic Jet Jt spectrum.
void SetDataPath(TString s)
Double_t Area() const
Definition: AliEmcalJet.h:123
TObjArray fClusterCollArray
cluster collection array
void setAssocPtBorders(int n, Double_t *c)
virtual AliVParticle * GetNextAcceptParticle()
Double_t GetRhoVal() const
void SetPeriodName(TString s)
double Double_t
Definition: External.C:58
const int debug
Definition: scanAll.C:15
void ExecOnce()
Perform steps needed to initialize the analysis.
AliJetContainer * GetJetContainer(Int_t i=0) const
void Terminate(Option_t *option)
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
AliClusterContainer * fCaloClustersCont
Tracks.
TH1 ** fHistClustersPt
Track eta spectrum.
Double_t fMinBinPt
min pt in histograms
TH1 **** fHistLogJTBg
Jt background.
JTJTEfficiency & operator=(const JTJTEfficiency &obj)
TCanvas * c
Definition: TestFitELoss.C:172
TProfile ** fhTrackingEfficiency
Track pT in jet.
Int_t fCentBin
!event centrality bin
TH1 *** fHistJetMulti
Jet eta distribution.
Int_t fNpttBins
Tracking efficiency.
TH1 *** fHistJetTracksPt
Multiplicity in jet.
TH1 **** fHistJTPtaNonInv
Associated pT vs. Jt in background.
Double_t getJt(AliVTrack *track, AliEmcalJet *jet, int reverse)
TH1 ** fHistTracksJt
Track pt spectrum.
Bool_t FillHistograms()
Function filling histograms.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Double_t getDiffR(double phi1, double phi2, double eta1, double eta2)
AliEmcalJet * GetLeadingJet(const char *opt="")
TProfile ** fHistJetsCorrPtVsNonCorr
Rho corrected Jet pt spectrum.
void SetRunNumber(Long64_t runnum)
TH1 *** fHistJetEta
Background pt distribution.
int Int_t
Definition: External.C:63
const AliVVertex * fPrimaryVertex
Event selection statistics.
Definition: External.C:204
AliParticleContainer * fTracksCont
Jets.
float Float_t
Definition: External.C:68
TDirectory * fEffDir[3]
void setCentBinBorders(int n, Double_t *c)
virtual AliVParticle * GetParticle(Int_t i=-1) const
TString fTrackArrayName
tracks array
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
void setTriggPtBorders(int n, Double_t *c)
Double_t fCent
!event centrality
TH1 *** fHistBgPt
Multiplicity in background cone.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Double_t CentBinBorders[10]
Vertex selection helper.
TH1 **** fHistLogJTPtaNonInv
Jet Jt spectrum.
TH1 *** fHistJetsPt
Leading jet pt spectrum.
TH1 ** fHistLeadingJetPt
Cluster pt spectrum.
AliVCluster * GetCluster(Int_t i) const
TH1 **** fHistJTPta_allNonInv
Logarithmic Jet Jt spectrum.
AliEmcalJet * GetNextAcceptJet()
TGraphErrors * fCorrection[20][20][20]
TObjArray fJetCollArray
jet collection array
TClonesArray * fTracks
Vertex found per event.
double GetCorrection(double pt, int icut, double cent) const
TH1 **** fHistJTBg
All particles Jt spectrum.
Double_t Pt() const
Definition: AliEmcalJet.h:102
TH1 *** fHistJetsCorrPt
Jet pt spectrum.
TH1 ** fHistTracksEta
Track jt spectrum.
TH1 **** fHistLogJTPta
Jet Jt spectrum.
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1 **** fHistLogJTBgNonInv
Jt background.
Double_t fVertex[3]
!event vertex
TH1 **** fHistJTBgNonInv
All particles Jt spectrum.
TH1 *** fHistBgMulti
Logarithmic Jt background.
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
AliAnalysisUtils * fVertexHelper
AliJ Efficiency.
Double_t Pz() const
Definition: AliEmcalJet.h:101
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
bool Bool_t
Definition: External.C:53
TH1 *** fHistBackgroundDone
Corrected versus raw jet pt.
AliVCluster * GetNextAcceptCluster()
Double_t M() const
Definition: AliEmcalJet.h:113
#define JUNUSED(expr)
TH1I * fHistEvtSelection
vertexZ inclusive
Definition: External.C:196
TProfile *** fHistPtaVsJt
Logarithmic Jt background.
Int_t fNbins
no. of pt bins
TH1 **** fHistJTPta
Background test.
TProfile *** fHistBgPtaVsJt
Associated pT vs. Jt in jet.