Coverity fixes
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetJTJT.cxx
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 "AliEmcalJet.h"
43 #include "AliRhoParameter.h"
44 #include "AliLog.h"
45 #include "AliJetContainer.h"
46 #include "AliParticleContainer.h"
47 #include "AliClusterContainer.h"
48 #include "AliPicoTrack.h"
49
50 #include "AliAnalysisTaskJetJTJT.h"
51
52
53 ClassImp(AliAnalysisTaskJetJTJT)
54
55         //________________________________________________________________________
56         AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT() : 
57                 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetJTJT", kTRUE),
58                 fHistTracksPt(0),
59                 fHistTracksJt(0),
60                 fHistClustersPt(0),
61                 fHistLeadingJetPt(0),
62                 fHistJetsPt(0),
63                 fHistBackgroundDone(0),
64                 fHistJTPta(0),
65                 fHistLogJTPta(0),
66                 fHistJTPta_all(0),
67                 fHistJTBg(0),
68                 fHistLogJTBg(0),
69                 fHistBgMulti(0),
70                 fHistBgPt(0),
71                 fHistJetEta(0),
72                 fHistJetMulti(0),
73                 fHistJetTracksPt(0),
74                 fhTrackingEfficiency(0),
75                 fNpttBins(1),
76                 fNptaBins(1),
77                 fJetsCont(0),
78                 //fJetsConts(0),
79                 //nJetsConts(0),
80                 fTracksCont(0),
81                 fCaloClustersCont(0),
82                 fTracks(0),
83                 fTrackArrayName("nonejk"),
84                 runPeriod(""),
85                 fEfficiency(0),
86                 debug(0)
87
88
89 {
90         // Default constructor.
91
92         fHistTracksPt       = new TH1*[fNcentBins];
93         fHistTracksJt       = new TH1*[fNcentBins];
94         fHistClustersPt     = new TH1*[fNcentBins];
95         fHistLeadingJetPt   = new TH1*[fNcentBins];
96         fHistJetsPt         = new TH1**[fNcentBins];
97         fHistBackgroundDone = new TH1**[fNcentBins];
98         fHistJTPta          = new TH1***[fNcentBins];
99         fHistLogJTPta       = new TH1***[fNcentBins];
100         fHistJTPta_all      = new TH1***[fNcentBins];
101         fHistJTBg           = new TH1***[fNcentBins];
102         fHistLogJTBg        = new TH1***[fNcentBins];
103         fHistBgMulti        = new TH1**[fNcentBins];
104         fHistBgPt           = new TH1**[fNcentBins];     
105         fHistJetEta         = new TH1**[fNcentBins];
106         fHistJetMulti       = new TH1**[fNcentBins];
107         fHistJetTracksPt     = new TH1**[fNcentBins];
108         fhTrackingEfficiency = new TProfile*[fNcentBins];
109         //CentBinBorders      = new Double_t[10];
110
111
112         for (Int_t i = 0; i < fNcentBins; i++) {
113                 fHistJTPta[i] = 0;
114                 fHistLogJTPta[i] = 0;
115                 fHistJTPta_all[i] = 0;
116                 fHistJTBg[i] = 0;
117                 fHistLogJTBg[i] = 0;
118                 fHistBackgroundDone[i] = 0;
119                 fHistTracksPt[i] = 0;
120                 fHistTracksJt[i] = 0;
121                 fHistClustersPt[i] = 0;
122                 fHistLeadingJetPt[i] = 0;
123                 fHistJetsPt[i] = 0;
124                 fHistBgMulti[i] = 0;
125                 fHistBgPt[i] = 0;  
126                 fHistJetEta[i] = 0; 
127                 fHistJetMulti[i] = 0;
128                 fHistJetTracksPt[i] = 0;
129                 fhTrackingEfficiency[i] = 0;
130         }
131
132         /*for(Int_t i = 0; i < nJetsConts; i++){
133           fJetsConts[i] = 0;
134           }*/
135         SetMakeGeneralHistograms(kTRUE);
136 }
137
138 //________________________________________________________________________
139 AliAnalysisTaskJetJTJT::AliAnalysisTaskJetJTJT(const char *name) : 
140         AliAnalysisTaskEmcalJet(name, kTRUE),
141         fHistTracksPt(0),
142         fHistTracksJt(0),
143         fHistClustersPt(0),
144         fHistLeadingJetPt(0),
145         fHistJetsPt(0),
146         fHistBackgroundDone(0),
147         fHistJTPta(0),
148         fHistLogJTPta(0),
149         fHistJTPta_all(0),
150         fHistJTBg(0),
151         fHistLogJTBg(0),
152         fHistBgMulti(0),
153         fHistBgPt(0),
154         fHistJetEta(0),
155         fHistJetMulti(0),
156         fHistJetTracksPt(0),
157         fhTrackingEfficiency(0),
158         fNpttBins(1),
159         fNptaBins(1),
160         fJetsCont(0),
161         //fJetsConts(0),
162         //nJetsConts(0),
163         fTracksCont(0),
164         fCaloClustersCont(0),
165         fTracks(0),
166         fTrackArrayName("nonejk"),
167         runPeriod(""),
168         fEfficiency(0),
169         debug(0)
170 {
171         // Standard constructor.
172         fHistTracksPt       = new TH1*[fNcentBins];
173         fHistTracksJt       = new TH1*[fNcentBins];
174         fHistClustersPt     = new TH1*[fNcentBins];
175         fHistLeadingJetPt   = new TH1*[fNcentBins];
176         fHistJetsPt         = new TH1**[fNcentBins];
177         fHistBackgroundDone = new TH1**[fNcentBins];    
178         fHistJTPta          = new TH1***[fNcentBins];   
179         fHistLogJTPta       = new TH1***[fNcentBins];   
180         fHistJTPta_all      = new TH1***[fNcentBins];   
181         fHistJTBg           = new TH1***[fNcentBins];   
182         fHistLogJTBg        = new TH1***[fNcentBins];   
183         fHistBgMulti        = new TH1**[fNcentBins];
184         fHistBgPt           = new TH1**[fNcentBins];
185         fHistJetEta         = new TH1**[fNcentBins];
186         fHistJetMulti       = new TH1**[fNcentBins];
187         fHistJetTracksPt     = new TH1**[fNcentBins];
188         fhTrackingEfficiency = new TProfile*[fNcentBins];
189         //CentBinBorders            = new Double_t[10];
190
191
192         for (Int_t i = 0; i < fNcentBins; i++) {
193                 fHistJTPta[i] = 0;
194                 fHistLogJTPta[i] = 0;
195                 fHistJTPta_all[i] = 0;
196                 fHistJTBg[i] = 0;
197                 fHistLogJTBg[i] = 0;
198                 fHistBackgroundDone[i] = 0;
199                 fHistTracksPt[i] = 0;
200                 fHistTracksJt[i] = 0;
201                 fHistClustersPt[i] = 0;
202                 fHistLeadingJetPt[i] = 0;
203                 fHistJetsPt[i] = 0;
204                 fHistBgMulti[i] = 0;        
205                 fHistBgPt[i] = 0;           
206                 fHistJetEta[i] = 0;         
207                 fHistJetMulti[i] = 0; 
208                 fHistJetTracksPt[i] = 0; 
209                 fhTrackingEfficiency[i] = 0;
210         }
211
212         /*for(Int_t i = 0; i < nJetsConts; i++){
213           fJetsConts[i] = 0;
214           }*/
215         SetMakeGeneralHistograms(kTRUE);
216 }
217
218 //________________________________________________________________________
219 AliAnalysisTaskJetJTJT::~AliAnalysisTaskJetJTJT()
220 {
221         // Destructor.
222 }
223
224
225 void AliAnalysisTaskJetJTJT::setCentBinBorders( int n, Double_t *c){
226         fNcentBins=n;  
227         if(debug > 0){
228                 cout << "AliAnalysisTaskJetJTJT::setCentBinBorders: " << endl;
229         }
230         for(int i= 0 ; i < fNcentBins; i++){
231                 CentBinBorders[i]= c[i];
232                 if(debug > 0)
233                         cout << CentBinBorders[i] << endl;
234         }       
235 }
236
237 void AliAnalysisTaskJetJTJT::setTriggPtBorders( int n, Double_t *c){
238         fNpttBins=n;  
239         if(debug > 0)
240                 cout << "AliAnalysisTaskJetJTJT::setTriggPtBorders: " << endl;
241         for(int i= 0 ; i < fNpttBins; i++){
242                 TriggPtBorders[i]= c[i];
243                 if(debug > 0)
244                         cout << TriggPtBorders[i] << endl;
245         }
246 }
247
248 void AliAnalysisTaskJetJTJT::setAssocPtBorders( int n, Double_t *c){
249         fNptaBins=n;  
250         if(debug > 0)
251                 cout << "AliAnalysisTaskJetJTJT::setAssocPtBorders: " << endl;
252         for(int i= 0 ; i < fNptaBins; i++){
253                 AssocPtBorders[i]= c[i];
254                 if(debug > 0)
255                         cout << AssocPtBorders[i] << endl;
256         }
257 }
258
259
260 //________________________________________________________________________
261 void AliAnalysisTaskJetJTJT::UserCreateOutputObjects()
262 {
263         // Create user output.
264
265         AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
266         if(debug > 0)
267                 cout << "Creating Histograms" << endl;
268
269         fJetsCont           = GetJetContainer(0);
270         /*for(int i = 0; i < nJetsConts ; i++){
271           fJetsConts[0]     = GetJetContainer(i);
272           }*/
273         /*if(fJetsCont) { //get particles and clusters connected to jets
274           fTracksCont       = fJetsCont->GetParticleContainer();
275           fCaloClustersCont = fJetsCont->GetClusterContainer();
276           } else {        //no jets, just analysis tracks and clusters
277           fTracksCont       = GetParticleContainer(0);
278           fCaloClustersCont = GetClusterContainer(0);
279           }*/
280         fTracksCont       = GetParticleContainer(0);
281         fCaloClustersCont = GetClusterContainer(0);
282         fTracksCont->SetClassName("AliVTrack");
283         fCaloClustersCont->SetClassName("AliAODCaloCluster");
284
285         TString histname;
286         //Int_t fMinBinJt = 0;
287         //Int_t fMaxBinJt = 5;
288
289         Int_t NBINSJt=150;
290         double LogBinsJt[NBINSJt+1], LimLJt=0.01, LimHJt=10;
291         double logBWJt = (TMath::Log(LimHJt)-TMath::Log(LimLJt))/(NBINSJt-1);
292         LogBinsJt[0] = 0;
293         for(int ij=1;ij<=NBINSJt;ij++) LogBinsJt[ij]=LimLJt*exp(ij*logBWJt);
294
295         int NBINSJtW=150;
296         double LimLJtW=TMath::Log(0.01), LimHJtW=TMath::Log(10);
297
298         //==== Efficiency ====
299         if(debug > 0)
300                 cout << "AliAnalysisTaskJetJTJT::UserCreateOutputObjects: Creating efficiency" << endl;
301         fEfficiency = new JTJTEfficiency;
302         fEfficiency->SetMode( 1 ); // 0:NoEff, 1:Period 2:RunNum 3:Auto
303         fEfficiency->SetDataPath("alien:///alice/cern.ch/user/d/djkim/legotrain/efficieny/data"); // Efficiency root file location local or alien
304
305
306
307         for (Int_t ic = 0; ic < fNcentBins; ic++) {
308                 if (fParticleCollArray.GetEntriesFast()>0) {
309                         histname = "fHistTracksPt_";
310                         histname += ic;
311                         fHistTracksPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt / 4);
312                         fHistTracksPt[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
313                         fHistTracksPt[ic]->GetYaxis()->SetTitle("counts");
314                         fOutput->Add(fHistTracksPt[ic]);
315                 }
316
317                 if (fParticleCollArray.GetEntriesFast()>0) {
318                         histname = "fHistTracksJt_";
319                         histname += ic;
320                         fHistTracksJt[ic] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
321                         fHistTracksJt[ic]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
322                         fHistTracksJt[ic]->GetYaxis()->SetTitle("counts");
323                         fOutput->Add(fHistTracksJt[ic]);
324                 }
325
326                 if (fParticleCollArray.GetEntriesFast()>0) {
327                         histname = "fhTrackingEfficiency_";
328                         histname += ic;
329                         fhTrackingEfficiency[ic] = new TProfile(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
330                         fhTrackingEfficiency[ic]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
331                         fhTrackingEfficiency[ic]->GetYaxis()->SetTitle("counts");
332                         fOutput->Add(fhTrackingEfficiency[ic]);
333                 }
334                 fHistJTPta[ic] = new TH1**[fNpttBins];
335                 fHistLogJTPta[ic] = new TH1**[fNpttBins];
336                 fHistJTPta_all[ic] = new TH1**[fNpttBins];
337                 fHistJTBg[ic] = new TH1**[fNpttBins];
338                 fHistLogJTBg[ic] = new TH1**[fNpttBins];
339                 fHistJetsPt[ic] = new TH1*[fNpttBins];
340                 fHistBackgroundDone[ic] = new TH1*[fNpttBins];
341                 fHistBgMulti[ic] = new TH1*[fNpttBins];
342                 fHistBgPt[ic] = new TH1*[fNpttBins];
343                 fHistJetEta[ic] = new TH1*[fNpttBins];
344                 fHistJetMulti[ic] = new TH1*[fNpttBins];
345                 fHistJetTracksPt[ic] = new TH1*[fNpttBins];
346                 for(Int_t j=0; j < fNpttBins; j++){
347                         fHistJTPta[ic][j] = new TH1*[fNptaBins];
348                         fHistLogJTPta[ic][j] = new TH1*[fNptaBins];
349                         fHistJTPta_all[ic][j] = new TH1*[fNptaBins];
350                         fHistJTBg[ic][j] = new TH1*[fNptaBins];
351                         fHistLogJTBg[ic][j] = new TH1*[fNptaBins];
352                         for(Int_t k=0; k < fNptaBins; k++){
353                                 fHistJTPta[ic][j][k] = 0;
354                                 fHistLogJTPta[ic][j][k] = 0;
355                                 fHistJTPta_all[ic][j][k] = 0;
356                                 fHistJTBg[ic][j][k] = 0;
357                                 fHistLogJTBg[ic][j][k] = 0;
358                         }
359                         fHistJetsPt[ic][j] = 0;
360                         fHistBackgroundDone[ic][j] = 0;
361                         fHistBgMulti[ic][j] = 0;
362                         fHistBgPt[ic][j] = 0;
363                         fHistJetEta[ic][j] = 0;
364                         fHistJetMulti[ic][j] =0;
365                         fHistJetTracksPt[ic][j] = 0;
366                 }
367
368
369                 if (fParticleCollArray.GetEntriesFast()>0) {
370                         for(Int_t iptt = 0 ; iptt <  fNpttBins; iptt++){
371                                 for(Int_t ipta = 0 ; ipta < fNptaBins; ipta++){
372                                         histname = "hJTPtaD00C";
373                                         //histname += ic;
374                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
375                                         if(debug > 1)
376                                                 cout << histname << endl;
377                                         fHistJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(),NBINSJt, LogBinsJt);
378                                         fHistJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
379                                         fHistJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
380                                         fOutput->Add(fHistJTPta[ic][iptt][ipta]);
381
382                                         histname = "hLogJTPtaD00C";
383                                         //histname += ic;
384                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
385                                         if(debug > 1)
386                                                 cout << histname << endl;
387                                         fHistLogJTPta[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
388                                         fHistLogJTPta[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
389                                         fHistLogJTPta[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
390                                         fOutput->Add(fHistLogJTPta[ic][iptt][ipta]);
391
392                                         histname = "hJTPta_allD00C";
393                                         //histname += ic;
394                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
395                                         if(debug > 1)
396                                                 cout << histname << endl;
397                                         fHistJTPta_all[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
398                                         fHistJTPta_all[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
399                                         fHistJTPta_all[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
400                                         fOutput->Add(fHistJTPta_all[ic][iptt][ipta]);
401
402                                         histname = "hJTBg";
403                                         //histname += ic;
404                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
405                                         if(debug > 1)
406                                                 cout << histname << endl;
407                                         fHistJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJt, LogBinsJt);
408                                         fHistJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("J_{T,track} (GeV/c)");
409                                         fHistJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
410                                         fOutput->Add(fHistJTBg[ic][iptt][ipta]);
411
412                                         histname = "hLogJTBg";
413                                         //histname += ic;
414                                         histname += Form("C%02dT%02dA%02d", ic, iptt, ipta);
415                                         if(debug > 1)
416                                                 cout << histname << endl;
417                                         fHistLogJTBg[ic][iptt][ipta] = new TH1F(histname.Data(), histname.Data(), NBINSJtW, LimLJtW, LimHJtW);
418                                         fHistLogJTBg[ic][iptt][ipta]->GetXaxis()->SetTitle("ln(J_{T,track}/ (GeV/c))");
419                                         fHistLogJTBg[ic][iptt][ipta]->GetYaxis()->SetTitle("counts");
420                                         fOutput->Add(fHistLogJTBg[ic][iptt][ipta]);
421
422
423                                 }
424                         }
425                 }
426
427                 if (fClusterCollArray.GetEntriesFast()>0) {
428                         histname = "fHistClustersPt_";
429                         histname += ic;
430                         fHistClustersPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
431                         fHistClustersPt[ic]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
432                         fHistClustersPt[ic]->GetYaxis()->SetTitle("counts");
433                         fOutput->Add(fHistClustersPt[ic]);
434                 }
435
436                 if (fJetCollArray.GetEntriesFast()>0) {
437                         histname = "fHistLeadingJetPt_";
438                         histname += ic;
439                         fHistLeadingJetPt[ic] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
440                         fHistLeadingJetPt[ic]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
441                         fHistLeadingJetPt[ic]->GetYaxis()->SetTitle("counts");
442                         fOutput->Add(fHistLeadingJetPt[ic]);
443
444                         for(Int_t iptt = 0 ; iptt <  fNpttBins; iptt++){
445
446                                 histname = "fHistJetsPt_";
447                                 histname += Form("C%02dT%02d", ic, iptt);
448                                 fHistJetsPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
449                                 fHistJetsPt[ic][iptt]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
450                                 fHistJetsPt[ic][iptt]->GetYaxis()->SetTitle("counts");
451                                 fOutput->Add(fHistJetsPt[ic][iptt]);
452
453                                 histname = "fHistBackgroundDone_";
454                                 histname += Form("C%02dT%02d", ic, iptt);;
455                                 fHistBackgroundDone[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 2, -1, 2);
456                                 fHistBackgroundDone[ic][iptt]->GetXaxis()->SetTitle("Number of jets");
457                                 fHistBackgroundDone[ic][iptt]->GetYaxis()->SetTitle("0 = not done, 1 = done");
458                                 fOutput->Add(fHistBackgroundDone[ic][iptt]);
459
460                                 histname = "fHistJetEta_";
461                                 histname += Form("C%02dT%02d", ic, iptt);
462                                 fHistJetEta[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, -5, 5);
463                                 fHistJetEta[ic][iptt]->GetXaxis()->SetTitle("#eta");
464                                 fHistJetEta[ic][iptt]->GetYaxis()->SetTitle("jets");
465                                 fOutput->Add(fHistJetEta[ic][iptt]);
466
467                                 histname = "fHistJetMulti_";
468                                 histname += Form("C%02dT%02d", ic, iptt);
469                                 fHistJetMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
470                                 fHistJetMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
471                                 fHistJetMulti[ic][iptt]->GetYaxis()->SetTitle("jets");
472                                 fOutput->Add(fHistJetMulti[ic][iptt]);
473
474                                 histname = "fHistBgMulti_";
475                                 histname += Form("C%02dT%02d", ic, iptt);
476                                 fHistBgMulti[ic][iptt] = new TH1F(histname.Data(), histname.Data(), 200, 0, 200);
477                                 fHistBgMulti[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
478                                 fHistBgMulti[ic][iptt]->GetYaxis()->SetTitle("Events");
479                                 fOutput->Add(fHistBgMulti[ic][iptt]);
480
481                                 histname = "fHistBgPt_";
482                                 histname += Form("C%02dT%02d", ic, iptt);
483                                 fHistBgPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins,fMinBinPt, fMaxBinPt/4);
484                                 fHistBgPt[ic][iptt]->GetXaxis()->SetTitle("Multiplicity");
485                                 fHistBgPt[ic][iptt]->GetYaxis()->SetTitle("jets");
486                                 fOutput->Add(fHistBgPt[ic][iptt]);
487
488                                 histname = "fHistJetTracksPt_";
489                                 //histname += ic;
490                                 histname += Form("C%02dT%02d", ic, iptt);
491                                 if(debug > 1)
492                                         cout << histname << endl;
493                                 fHistJetTracksPt[ic][iptt] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt/4);
494                                 fHistJetTracksPt[ic][iptt]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
495                                 fHistJetTracksPt[ic][iptt]->GetYaxis()->SetTitle("counts");
496                                 fOutput->Add(fHistJetTracksPt[ic][iptt]);
497                         }
498
499                         /*
500                            if (!(GetJetContainer()->GetRhoName().IsNull())) {
501                            histname = "fHistJetsCorrPtArea_";
502                            histname += i;
503                            fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
504                            fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
505                            fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
506                            fOutput->Add(fHistJetsCorrPtArea[i]);
507                            }
508                            */
509                 }
510         }
511
512         PostData(1, fOutput); // Post data for ALL output slots > 0 here.
513 }
514
515 //________________________________________________________________________
516 Bool_t AliAnalysisTaskJetJTJT::FillHistograms()
517 {
518         // Fill histograms.
519
520         AliCentrality *aliCent = InputEvent()->GetCentrality();
521         fCentBin = 0;
522         if (aliCent) {
523                 //fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
524                 fCent = aliCent->GetCentralityPercentile("V0M");
525                 /*if(debug > 0){
526                   cout << "Centrality " << fCent << endl;
527                   }*/
528                 for(int ic = 0; ic < fNcentBins; ic++){
529                         /*if(debug > 0){
530                           cout << "ic: " << ic << " / " << fNcentBins << endl;
531                           cout << "Centrality bin " << fCentBin << endl;
532                           cout << "Border: " << CentBinBorders[ic] << endl;
533                           } */
534                         if(fCent > CentBinBorders[ic]){
535                                 fCentBin = ic;
536                         }
537                 }
538                 //cout << "Centrality bin: " << fCentBin << endl;
539         } else {
540                 AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
541                 fCentBin = 3;
542         }
543         int fHadronSelectionCut = 5; //5=Hybrid cut
544         if (fTracksCont) {
545                 AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
546                 while(track) {
547                         double ptt = track->Pt();
548
549                         //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
550                         //AliJBaseTrack *triggTr = (AliJBaseTrack*)fchargedHadronList->At(i);
551
552                         double effCorr = 1./fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
553                         //double effCorr = 1.;
554                         fhTrackingEfficiency[fCentBin]->Fill( ptt, 1./effCorr );
555                         //triggTr->SetTrackEff( 1./effCorr );
556                         //<<<<<<<<<<<< Efficiency >>>>>>>>>>>
557
558                         if(ptt > 0 && 1.0/ptt > 0){
559                                 fHistTracksPt[fCentBin]->Fill(ptt,1./ptt*effCorr); 
560                         }
561
562
563                         track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
564                 }
565         }
566
567         if (fCaloClustersCont) {
568                 AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
569                 while(cluster) {
570                         TLorentzVector nPart;
571                         cluster->GetMomentum(nPart, fVertex);
572                         fHistClustersPt[fCentBin]->Fill(nPart.Pt());
573
574                         cluster = fCaloClustersCont->GetNextAcceptCluster();
575                 }
576         }
577
578         Int_t fPttBin, fPtaBin;
579         fPtaBin = 0;
580
581
582         if (fJetsCont) {
583                 //Int_t Njets = fJetsCont->GetNJets();
584                 Int_t Njets = 0;
585                 if(debug > 1){
586                         cout << "Number of Jets: " << Njets << endl;
587                 }
588
589                 //Make a array to hold jets to be tested in background jT
590                 Float_t jetPhis[200] = {};
591                 Float_t jetEtas[200] = {};
592                 AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
593                 Int_t ij = 0;
594                 while(jet) {
595                         //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
596                         if(jet->Pt() > 5){    //Only consider jets with pT > 5 GeV
597                                 jetPhis[ij] = jet->Phi();
598                                 jetEtas[ij] = jet->Eta();
599                                 ij++;
600                                 Njets++;
601                                 if(debug > 1)
602                                         cout << "i: " << ij << " jetPhi: " << jetPhis[ij] << " jetEta: " << jetEtas[ij] << endl;
603                         }else{
604                                 //jetPhis[ij] = 100;
605                                 //jetEtas[ij] = 100;
606                                 //Njets--;
607                                 if(debug > 1)
608                                         cout << "jetPt: " << jet->Pt() << " jetPhi: " << jet->Phi() << " jetEta: " << jet->Eta() << endl;
609                         }
610                         //i++; 
611                         jet = fJetsCont->GetNextAcceptJet();
612                 }
613
614
615                 //fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
616                 jet = fJetsCont->GetNextAcceptJet(0); 
617                 while(jet) {
618                         if(jet->Pt() > 5){
619                                 if(jet->Eta() < -0.4 || jet->Eta() > 0.4){
620                                         if(debug > 0)
621                                                 cout << "Jet outside eta range, Eta: " << jet->Eta() << endl;
622                                         jet = fJetsCont->GetNextAcceptJet();
623                                         continue;
624                                 }
625                                 //cout << "Jet found " << ij << " pt: " << jet->Pt() << endl;
626                                 //Get the trigger pT bin
627                                 fPttBin = 0;
628                                 for(int iptt = 0 ; iptt < fNpttBins; iptt++){
629                                         if(jet->Pt() > TriggPtBorders[iptt]){
630                                                 fPttBin = iptt;
631                                         }
632
633                                 }
634                                 fHistJetEta[fCentBin][fPttBin]->Fill(jet->Eta());
635                                 if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
636                                         fHistJetsPt[fCentBin][fPttBin]->Fill(jet->Pt(),1.0/jet->Pt());  //Fill jet dN/(pT dpT)
637                                 }
638
639                                 /*if (fHistJetsCorrPtArea[fCentBin]) {
640                                   Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
641                                   fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
642                                   }*/
643                                 //Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
644
645
646                                 Int_t nTrack = jet->GetNumberOfTracks();
647                                 if (debug > 0)                  
648                                         cout << "Number of tracks " << nTrack << " Jet Pt: " << jet->Pt() << endl;
649                                 fHistJetMulti[fCentBin][fPttBin]->Fill(nTrack);
650                                 for(Int_t it = 0; it < nTrack; it++ ){
651                                         AliVParticle *track = (AliVParticle*)jet->TrackAt( it, fTracks );
652                                         if( !track ){
653                                                 cout << "No Track found" << endl;
654                                                 continue;
655                                         }
656                                         fPtaBin = 0; //Get the associated pT bin
657                                         for(int ipta = 0 ; ipta < fNptaBins; ipta++){
658                                                 if(track->Pt() > AssocPtBorders[ipta]){
659                                                         fPtaBin = ipta;
660                                                 }
661                                         }
662                                         fHistJetTracksPt[fCentBin][fPttBin]->Fill(track->Pt());
663                                         if(debug > 2)
664                                                 cout << "Filling fHistJetTracksPt C" << fCentBin << " T" << fPttBin << endl;
665                                         /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz(); // p_track dot p_jet
666                                         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz()); // track pT norm
667                                         Float_t normproduct = constp*jetp; // jet pT norm times track pT norm
668                                         Float_t costheta2 = dotproduct/normproduct; 
669                                         //Float_t sintheta = sqrt(1-costheta2*costheta2);
670                                         Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
671                                         Float_t jt = getJt(track,jet,0);
672                                         double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
673                                         //double effCorr = 1.;
674                                         if(jt > 0 && 1.0/jt > 0){
675                                                 fHistTracksJt[fCentBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
676                                                 fHistJTPta[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr); //Fill dN/(djT jT)
677                                                 fHistLogJTPta[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr); //Fill logarithmic dN/(djT jT)
678                                         }
679                                         if(debug > 1)
680                                                 cout << "Filling JT C" << fCentBin << "T" <<  fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
681                                 }
682
683                                 //Get Jet azimuth and rapidity of jet
684                                 Float_t jetAngle = jet->Phi();
685                                 Float_t jetRap = jet->Eta();
686
687                                 //Rotate jet angle for background cone
688                                 Float_t rotatedAngle = jetAngle+TMath::Pi()/2;
689                                 if(rotatedAngle > TMath::Pi()*2){
690                                         rotatedAngle = rotatedAngle- TMath::Pi()*2;
691                                 }
692                                 Float_t jetArea = jet->Area();
693                                 Float_t testRadius = TMath::Sqrt(jetArea/TMath::Pi());
694
695                                 Bool_t doBg = 1;
696
697                                 //Test if there are jets in the background test cone
698                                 for(int i_j = 0; i_j < Njets; i_j++){
699                                         Float_t diffR = TMath::Sqrt(TMath::Power(jetPhis[i_j]-rotatedAngle,2)+TMath::Power(jetEtas[i_j]-jetRap,2));
700                                         if(debug > 1){
701                                                 cout << "i_j: " << i_j << " JetPhi: " << jetPhis[i_j] << " jetEta: " << jetEtas[i_j] << endl;
702                                                 cout << "DiffR: " << diffR << " doBG: " << doBg <<endl;
703                                         }
704                                         if(diffR < testRadius *2){ //Jets muts be at least 2*cone radius away from the background cone axis
705                                                 doBg =0;
706                                                 break;
707                                         }
708
709                                 }
710
711                                 // Do jT for all particles in respect to jet axis
712                                 if (fTracksCont) {
713                                         int counter = 0;
714                                         AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
715                                         while(track) {
716                                                 /*Float_t dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
717                                                   Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
718                                                   Float_t normproduct = constp*jetp;
719                                                   Float_t costheta2 = dotproduct/normproduct;
720                                                 //Float_t sintheta = sqrt(1-costheta2*costheta2);
721                                                 Float_t jt = constp*sqrt(1-costheta2*costheta2);*/
722                                                 Double_t jt = getJt(track,jet,0);
723                                                 double effCorr = 1./fEfficiency->GetCorrection(track->Pt(), fHadronSelectionCut, fCent);  // here you generate warning if ptt>30
724                                                 if(jt > 0 && 1.0/jt > 0){
725                                                         fHistJTPta_all[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
726                                                 }
727                                                 for(int ipta = 0 ; ipta < fNptaBins; ipta++){
728                                                         if(track->Pt() > AssocPtBorders[ipta]){
729                                                                 fPtaBin = ipta;
730                                                         }
731                                                 }
732
733                                                 //If background is to be filled
734                                                 if(doBg){
735                                                         Float_t diffR = TMath::Sqrt(TMath::Power(track->Phi()-rotatedAngle,2)+TMath::Power(track->Eta()-jetRap,2));
736                                                         //Particles in the rotated cone
737                                                         if(diffR < testRadius){
738                                                                 counter++;
739                                                                 fHistBgPt[fCentBin][fPttBin]->Fill(track->Pt());
740                                                                 /*dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
741                                                                   constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
742                                                                   normproduct = constp*jetp;
743                                                                   costheta2 = dotproduct/normproduct;
744                                                                 //sintheta = sqrt(1-costheta2*costheta2);
745                                                                 jt = constp*sqrt(1-costheta2*costheta2);*/
746                                                                 jt = getJt(track,jet,1);
747                                                                 if(jt > 0 && 1.0/jt > 0){
748                                                                         fHistJTBg[fCentBin][fPttBin][fPtaBin]->Fill(jt,1.0/jt*effCorr);
749                                                                         fHistLogJTBg[fCentBin][fPttBin][fPtaBin]->Fill(TMath::Log(jt),1.0/jt*effCorr);
750                                                                 }
751                                                                 if(debug > 1)
752                                                                         cout << "Filling Background C" << fCentBin << "T" <<  fPttBin << "A" << fPtaBin << " jt:" << jt << " with " << 1.0/jt*effCorr << endl;
753                                                                 //Fill background jT
754
755                                                         }
756                                                 }
757                                                 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
758                                         }
759                                         if(doBg){
760                                                 fHistBgMulti[fCentBin][fPttBin]->Fill(counter);
761                                         }
762                                 }
763                                 if(doBg){
764                                         fHistBackgroundDone[fCentBin][fPttBin]->Fill(1);
765                                 }else{
766                                         fHistBackgroundDone[fCentBin][fPttBin]->Fill(0);
767                                 }
768
769                         }
770                         jet = fJetsCont->GetNextAcceptJet(); 
771
772                 }
773                 jet = fJetsCont->GetLeadingJet();
774                 if(jet){
775                         if(jet->Pt() > 0 && 1.0/jet->Pt() > 0){
776                                 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt(),1.0/jet->Pt());
777                         }
778                 }
779         }
780
781         CheckClusTrackMatching();
782
783         return kTRUE;
784 }
785
786
787 //-----------------------------------------------------------------------
788 Double_t AliAnalysisTaskJetJTJT::getJt(AliVTrack *track, AliEmcalJet *jet,int reverse){
789         Float_t dotproduct = 0;
790         Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
791         if(reverse){
792                 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
793         } else{
794                 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
795         }
796         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
797         Float_t normproduct = constp*jetp;
798         Float_t costheta2 = dotproduct/normproduct;
799         //Float_t sintheta = sqrt(1-costheta2*costheta2);
800         Float_t jt = constp*sqrt(1-costheta2*costheta2);
801         return jt;
802 }
803
804 Double_t AliAnalysisTaskJetJTJT::getJt(AliVParticle *track, AliEmcalJet *jet,int reverse){
805         Float_t dotproduct = 0;
806         Float_t jetp = sqrt(jet->Px()*jet->Px()+jet->Py()*jet->Py()+jet->Pz()*jet->Pz()); //Jet pT norm
807         if(reverse){
808                 dotproduct = -track->Px()*jet->Py()+track->Py()*jet->Px()+track->Pz()*jet->Pz();
809         } else{
810                 dotproduct = track->Px()*jet->Px()+track->Py()*jet->Py()+track->Pz()*jet->Pz();
811         }
812         Float_t constp = sqrt(track->Px()*track->Px()+track->Py()*track->Py()+track->Pz()*track->Pz());
813         Float_t normproduct = constp*jetp;
814         Float_t costheta2 = dotproduct/normproduct;
815         //Float_t sintheta = sqrt(1-costheta2*costheta2);
816         Float_t jt = constp*sqrt(1-costheta2*costheta2);
817         return jt;
818 }
819
820 //________________________________________________________________________
821 void AliAnalysisTaskJetJTJT::CheckClusTrackMatching()
822 {
823
824         if(!fTracksCont || !fCaloClustersCont)
825                 return;
826
827         Double_t deta = 999;
828         Double_t dphi = 999;
829
830         //Get closest cluster to track
831         AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
832         while(track) {
833                 //Get matched cluster
834                 Int_t emc1 = track->GetEMCALcluster();
835                 if(fCaloClustersCont && emc1>=0) {
836                         AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
837                         if(clusMatch) {
838                                 AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
839                                 //fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
840                         }
841                 }
842                 track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
843         }
844
845         //Get closest track to cluster
846         AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
847         while(cluster) {
848                 TLorentzVector nPart;
849                 cluster->GetMomentum(nPart, fVertex);
850                 fHistClustersPt[fCentBin]->Fill(nPart.Pt());
851
852                 //Get matched track
853                 AliVTrack *mt = NULL;      
854                 AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
855                 if(acl) {
856                         if(acl->GetNTracksMatched()>1)
857                                 mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
858                 }
859                 else {
860                         AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
861                         Int_t im = ecl->GetTrackMatchedIndex();
862                         if(fTracksCont && im>=0) {
863                                 mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
864                         }
865                 }
866                 if(mt) {
867                         AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
868                         //fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
869
870                         //debugging
871                         /*
872                            if(mt->IsEMCAL()) {
873                            Int_t emc1 = mt->GetEMCALcluster();
874                            Printf("current id: %d  emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
875                            AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
876                            AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
877                            Printf("deta: %f dphi: %f",deta,dphi);
878                            }
879                            */    
880                 }
881                 cluster = fCaloClustersCont->GetNextAcceptCluster();
882         }
883 }
884
885 //________________________________________________________________________
886 /*AliJetContainer* AliAnalysisTaskJetJTJT::AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius) {
887
888   AliAnalysisTaskEmcalJet::ExecOnce();
889   nJetsConts++;
890   AliJetContainer *cont = 0x0;
891   cont = AliAnalysisTaskEmcalJet::AddJetContainer(n,defaultCutType,jetRadius);
892   return cont;
893   }*/
894
895
896 //________________________________________________________________________
897 void AliAnalysisTaskJetJTJT::ExecOnce() {
898
899         if(debug > 0){
900         cout << "AliAnalysisTaskJetJTJT::ExecOnce(): " << endl;
901         cout << "Get fTracks from " << fTrackArrayName.Data() << endl;
902         }
903         fTracks =dynamic_cast <TClonesArray*>(InputEvent()->FindListObject( fTrackArrayName.Data() ));
904
905         AliAnalysisTaskEmcalJet::ExecOnce();
906         if(debug > 1)
907                 cout << "Efficiency: Set Run Period Name " << runPeriod << endl;
908         fEfficiency->SetPeriodName(runPeriod);
909         if(debug > 1)
910                 cout << "Efficiency: Set Run number " << InputEvent()->GetRunNumber() << endl;
911         fEfficiency->SetRunNumber( InputEvent()->GetRunNumber() ); //TODO Get run Number
912         if(debug > 1)
913                 cout << "Efficiency: Load()" << endl;
914         fEfficiency->Load();
915         if(debug > 1)
916                 cout << "fEfficiency loaded" << endl;
917
918
919
920         if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
921         if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
922         if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
923
924 }
925
926 //________________________________________________________________________
927 Bool_t AliAnalysisTaskJetJTJT::Run()
928 {
929         // Run analysis code here, if needed. It will be executed before FillHistograms().
930
931         return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
932 }
933
934 //________________________________________________________________________
935 void AliAnalysisTaskJetJTJT::Terminate(Option_t *) 
936 {
937         // Called once at the end of the analysis.
938 }
939
940
941
942 //________________________________________________________________________
943 JTJTEfficiency::JTJTEfficiency():
944         fMode(kAuto),
945         fPeriod(-1),
946         fDataPath(""),
947         fName(""),
948         fPeriodStr(""),
949         fMCPeriodStr(""),
950         fRunNumber(0),
951         fTag(""),
952         fInputRootName(""),
953         fInputRoot(NULL),
954         fCentBinAxis(0x0)
955 {
956   for (Int_t i=0; i < 3; i++) fEffDir[i] = 0;
957
958 }
959
960 JTJTEfficiency::JTJTEfficiency(const JTJTEfficiency& obj) :
961   fMode(obj.fMode),
962   fPeriod(obj.fPeriod),
963   fDataPath(obj.fDataPath),
964   fName(obj.fName),
965   fPeriodStr(obj.fPeriodStr),
966   fMCPeriodStr(obj.fMCPeriodStr),
967   fRunNumber(obj.fRunNumber),
968   fTag(obj.fTag),
969   fInputRootName(obj.fInputRootName),
970   fInputRoot(obj.fInputRoot),
971   fCentBinAxis(obj.fCentBinAxis)
972 {
973         // copy constructor TODO: handling of pointer members
974         JUNUSED(obj);
975   for (Int_t i=0; i < 3; i++) fEffDir[i] = obj.fEffDir[i];
976 }
977
978 JTJTEfficiency& JTJTEfficiency::operator=(const JTJTEfficiency& obj){
979         // equal sign operator TODO: content
980         JUNUSED(obj);
981         return *this;
982 }
983
984
985 //________________________________________________________________________
986 double JTJTEfficiency::GetCorrection( double pt, int icut , double cent ) const {
987         if( fMode == kNotUse ) return 1;
988         int icent = fCentBinAxis->FindBin( cent ) -1 ;
989         if( icent < 0 || icent > fCentBinAxis->GetNbins()-1 ) {
990                 cout<<"J_WARNING : Centrality "<<cent<<" is out of CentBinBorder"<<endl;
991                 return 1;
992         }
993         // TODO error check for icent;
994         int ivtx = 0;
995         if( ! fCorrection[ivtx][icent][icut] ) {
996                 cout<<"J_WARNING : No Eff Info "<<pt<<"\t"<<icut<<"\t"<<cent<<"\t"<<icent<<endl;
997                 return 1;
998         }
999         TGraphErrors * gr = fCorrection[ivtx][icent][icut];
1000         //=== TEMPERORY SETTING. IT will be removed soon.
1001         if( pt > 30 ) pt = 30; // Getting eff of 30GeV for lager pt
1002         double cor = gr->Eval(pt);
1003         if ( cor < 0.2 ) cor = 0.2;
1004         return cor;
1005 }
1006
1007
1008 TString JTJTEfficiency::GetEffName() {
1009         /*
1010            1. kNotUse : no Load, efficiency is 1 always
1011            2. has fInputRootName : Load that or crash
1012            3. has fName : Load fName [+runnumber] or crash
1013            4. has runnumber : Find Good MC period from AliJRunTable, or crash
1014            3. has period : Find Good MC period from AliJRunTable, or crash
1015
1016            }
1017            */
1018 if(fPeriodStr == "LHC10b"){
1019         fInputRootName = "Eff--LHC10b-LHC10d1-0-.root";
1020 }
1021 if(fPeriodStr == "LHC10c"){
1022         fInputRootName = "Eff--LHC10c-LHC10d4-0-.root";
1023 }
1024 if(fPeriodStr == "LHC10d"){
1025         fInputRootName = "Eff--LHC10d-LHC10f6a-0-.root";
1026 }
1027 if(fPeriodStr == "LHC10e"){
1028         fInputRootName = "Eff--LHC10e-LHC10e20-0-.root";
1029 }
1030 if(fPeriodStr == "LHC10h"){
1031         fInputRootName = "Eff--LHC10h-LHC11a10a_bis-0-.root";
1032 }
1033 if(fPeriodStr == "LHC11a"){
1034         fInputRootName = "Eff--LHC11a-LHC11b10a-0-.root";
1035 }
1036 if(fPeriodStr == "LHC13b"){
1037         fInputRootName = "Eff--LHC13b-LHC13b2-efix_p1-0-.root";
1038 }
1039
1040 if(fPeriodStr == "LHC13c"){
1041         fInputRootName = "Eff--LHC13c-LHC13b2-efix_p1-0-.root";
1042 }
1043 if(fPeriodStr == "LHC13d"){
1044         fInputRootName = "Eff--LHC13d-LHC13b2-efix_p1-0-.root";
1045 }
1046 if(fPeriodStr == "LHC13e"){
1047         fInputRootName = "Eff--LHC13e-LHC13b2-efix_p1-0-.root";
1048 }
1049
1050 return fInputRootName;
1051 }
1052
1053 TString JTJTEfficiency::GetEffFullName() {
1054         GetEffName();
1055         fInputRootName = fDataPath + "/" + fInputRootName;
1056         return fInputRootName;
1057 }
1058
1059
1060 //________________________________________________________________________
1061 bool JTJTEfficiency::Load(){
1062         // Load Efficiency File based on fMode
1063         if( fMode == kNotUse ) {
1064                 cout<<"J_WARNING : Eff Mode is \"NOTUSE\". eff is 1 !!!"<<endl;
1065                 return true;
1066         }
1067         GetEffFullName();
1068         if (TString(fInputRootName).BeginsWith("alien:"))  TGrid::Connect("alien:");
1069         fInputRoot = TFile::Open( fInputRootName);
1070         //fInputRoot = new TFile( fInputRootName,"READ");
1071         if( !fInputRoot ) {
1072                 cout << "J_ERROR : %s does not exist" << fInputRootName << endl;
1073                 return false;
1074         }
1075
1076         //fEffDir[0] = (TDirectory*)fInputRoot->Get("EffRE");
1077         ///fEffDir[1] = (TDirectory*)fInputRoot->Get("EffMC");
1078         fEffDir[2] = (TDirectory*)fInputRoot->Get("Efficiency");
1079         //iif( fEffDir[0] && fEffDir[1] && fEffDir[2] )
1080         if( !fEffDir[2] )
1081         {
1082                 cout << "J_ERROR : Directory EFF is not exist"<<endl;
1083                 return false;
1084         }
1085
1086         fCentBinAxis = (TAxis*)fEffDir[2]->Get("CentralityBin");
1087         if( !fCentBinAxis ){
1088                 cout << "J_ERROR : No CentralityBin in directory" << endl;
1089                 gSystem->Exit(1);
1090         }
1091
1092
1093         int nVtx = 1;
1094         int nCentBin = fCentBinAxis->GetNbins();
1095         for( int ivtx=0;ivtx<nVtx;ivtx++ ){
1096                 for( int icent=0;icent<nCentBin;icent++ ){
1097                         for( int icut=0;icut<kJNTrackCuts;icut++ ){
1098                                 fCorrection[ivtx][icent][icut]
1099                                         = (TGraphErrors*) fEffDir[2]->Get(Form("gCor%02d%02d%02d", ivtx,icent,icut));
1100                                 //cout<<"J_LOG : Eff graph - "<<Form("gCor%02d%02d%02d", ivtx,icent,icut)<<" - "<<g<<endl;
1101                         }
1102                 }
1103         }
1104         cout<<"J_LOG : Eff file is "<<fInputRootName<<endl;
1105         cout<<"J_LOG : Eff Cent Bins are ";
1106         for( int i=0;i<=nCentBin;i++ ){
1107                 //cout<<fCentBinAxis->GetXbins()->At(i)<<" ";
1108         }
1109         //cout<<endl;
1110         return true;
1111 }
1112
1113
1114