coverity fix
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
1
2 // ******************************************
3 // This task computes several jet observables like 
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and 
6 // correlation strength distribution of particles inside jets.    
7 // Author: lcunquei@cern.ch
8 // *******************************************
9
10
11 /**************************************************************************
12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13  *                                                                        *
14  * Author: The ALICE Off-line Project.                                    *
15  * Contributors are mentioned in the code where appropriate.              *
16  *                                                                        *
17  * Permission to use, copy, modify and distribute this software and its   *
18  * documentation strictly for non-commercial purposes is hereby granted   *
19  * without fee, provided that the above copyright notice appears in all   *
20  * copies and that both the copyright notice and this permission notice   *
21  * appear in the supporting documentation. The authors make no claims     *
22  * about the suitability of this software for any purpose. It is          *
23  * provided "as is" without express or implied warranty.                  *
24  **************************************************************************/
25
26
27 #include "TChain.h"
28 #include "TTree.h"
29 #include "TMath.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH3F.h"
33 #include "THnSparse.h"
34 #include "TCanvas.h"
35
36 #include "AliLog.h"
37
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliVEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliCentrality.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliInputEventHandler.h"
47 #include "AliAODJetEventBackground.h"
48 #include "AliAODMCParticle.h"
49 #include "AliAnalysisTaskFastEmbedding.h"
50 #include "AliAODEvent.h"
51 #include "AliAODHandler.h"
52 #include "AliAODJet.h"
53
54 #include "AliAnalysisTaskJetCore.h"
55
56 using std::cout;
57 using std::endl;
58
59 ClassImp(AliAnalysisTaskJetCore)
60
61 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
62 AliAnalysisTaskSE(),
63 fESD(0x0),
64 fAODIn(0x0),
65 fAODOut(0x0),
66 fAODExtension(0x0),
67 fBackgroundBranch(""),
68 fNonStdFile(""),
69 fIsPbPb(kTRUE),
70 fOfflineTrgMask(AliVEvent::kAny),
71 fMinContribVtx(1),
72 fVtxZMin(-10.),
73 fVtxZMax(10.),
74 fEvtClassMin(0),
75 fEvtClassMax(4),
76 fFilterMask(0),
77 fFilterMaskBestPt(0),
78 fFilterType(0),
79 fRadioFrac(0.2),
80 fMinDist(0.1),
81 fCentMin(0.),
82 fCentMax(100.),
83 fNInputTracksMin(0),
84 fNInputTracksMax(-1),
85 fAngStructCloseTracks(0),
86 fCheckMethods(0),
87 fDoEventMixing(0), 
88 fFlagPhiBkg(0),
89 fFlagEtaBkg(0),
90 fFlagJetHadron(0),
91 fFlagRandom(0),
92 fFlagOnlyRecoil(0),
93 fFlagOnlyHardest(1),
94 fTrackTypeRec(kTrackUndef),
95 fRPAngle(0),
96 fNRPBins(50),
97 fSemigoodCorrect(0),
98 fHolePos(4.71),
99 fHoleWidth(0.2),
100 fJetEtaMin(-.5),
101 fJetEtaMax(.5),
102 fNevents(0),
103 fTindex(0),
104 fTrigBufferIndex(0),
105 fCountAgain(0), 
106 fJetPtMin(20.),
107 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
108 fJetPtFractionMin(0.5),
109 fNMatchJets(4),
110 fMatchMaxDist(0.8),
111 fKeepJets(kFALSE),
112 fRunAnaAzimuthalCorrelation(kFALSE),
113 fkNbranches(2),
114 fkEvtClasses(12),
115 fOutputList(0x0),
116 fbEvent(kTRUE),
117 fHistEvtSelection(0x0),
118 fhnDeltaR(0x0),
119 fhnMixedEvents(0x0),
120 fh2JetCoreMethod1C10(0x0),
121 fh2JetCoreMethod2C10(0x0),
122 fh2JetCoreMethod1C20(0x0),
123 fh2JetCoreMethod2C20(0x0),
124 fh2JetCoreMethod1C30(0x0),
125 fh2JetCoreMethod2C30(0x0),
126 fh2JetCoreMethod1C60(0x0),
127 fh2JetCoreMethod2C60(0x0),
128 fh3JetTrackC3060(0x0),
129 fh3JetTrackC20(0x0),
130 fh2AngStructpt1C10(0x0),
131 fh2AngStructpt2C10(0x0),
132 fh2AngStructpt3C10(0x0),
133 fh2AngStructpt4C10(0x0),
134 fh2AngStructpt1C20(0x0),
135 fh2AngStructpt2C20(0x0),
136 fh2AngStructpt3C20(0x0),
137 fh2AngStructpt4C20(0x0),    
138 fh2AngStructpt1C30(0x0),
139 fh2AngStructpt2C30(0x0),
140 fh2AngStructpt3C30(0x0),
141 fh2AngStructpt4C30(0x0),   
142 fh2AngStructpt1C60(0x0),
143 fh2AngStructpt2C60(0x0),
144 fh2AngStructpt3C60(0x0),
145 fh2AngStructpt4C60(0x0),
146 fh2Ntriggers(0x0),
147 fh2Ntriggers2C10(0x0),
148 fh2Ntriggers2C20(0x0), 
149 fh3JetDensity(0x0),
150 fh3JetDensityA4(0x0),
151 fh2RPJetsC10(0x0),
152 fh2RPJetsC20(0x0),
153 fh2RPTC10(0x0),
154 fh2RPTC20(0x0), 
155 fHJetSpec(0x0),
156 fhTTPt(0x0),
157 fHJetPhiCorr(0x0)
158
159  
160 {
161    // default Constructor
162
163
164  // Trigger buffer.
165    for(Int_t i=0; i<10; i++) {
166                 for(Int_t j=0; j<6; j++) {
167                         fTrigBuffer[i][j]=0;
168                 }                               
169         }       
170
171
172
173
174
175    fJetBranchName[0] = "";
176    fJetBranchName[1] = "";
177
178    fListJets[0] = new TList;
179    fListJets[1] = new TList;
180 }
181
182 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
183 AliAnalysisTaskSE(name),
184 fESD(0x0),
185 fAODIn(0x0),
186 fAODOut(0x0),
187 fAODExtension(0x0),
188 fBackgroundBranch(""),
189 fNonStdFile(""),
190 fIsPbPb(kTRUE),
191 fOfflineTrgMask(AliVEvent::kAny),
192 fMinContribVtx(1),
193 fVtxZMin(-10.),
194 fVtxZMax(10.),
195 fEvtClassMin(0),
196 fEvtClassMax(4),
197 fFilterMask(0),
198 fFilterMaskBestPt(0),
199 fFilterType(0),
200 fRadioFrac(0.2),
201 fMinDist(0.1),
202 fCentMin(0.),
203 fCentMax(100.),
204 fNInputTracksMin(0),
205 fNInputTracksMax(-1),
206 fAngStructCloseTracks(0),
207 fCheckMethods(0),
208 fDoEventMixing(0),
209 fFlagPhiBkg(0),
210 fFlagEtaBkg(0),
211 fFlagJetHadron(0),
212 fFlagRandom(0),
213 fFlagOnlyRecoil(0),
214 fFlagOnlyHardest(1),
215 fTrackTypeRec(kTrackUndef),
216 fRPAngle(0),
217 fNRPBins(50),
218 fSemigoodCorrect(0),
219 fHolePos(4.71),
220 fHoleWidth(0.2),
221 fJetEtaMin(-.5),
222 fJetEtaMax(.5),
223 fNevents(0),
224 fTindex(0),
225 fTrigBufferIndex(0),
226 fCountAgain(0),
227 fJetPtMin(20.),
228 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
229 fJetPtFractionMin(0.5),
230 fNMatchJets(4),
231 fMatchMaxDist(0.8),
232 fKeepJets(kFALSE),
233 fRunAnaAzimuthalCorrelation(kFALSE),
234 fkNbranches(2),
235 fkEvtClasses(12),
236 fOutputList(0x0),
237 fbEvent(kTRUE),
238 fHistEvtSelection(0x0),
239 fhnDeltaR(0x0),
240 fhnMixedEvents(0x0),
241 fh2JetCoreMethod1C10(0x0),
242 fh2JetCoreMethod2C10(0x0),
243 fh2JetCoreMethod1C20(0x0),
244 fh2JetCoreMethod2C20(0x0),
245 fh2JetCoreMethod1C30(0x0),
246 fh2JetCoreMethod2C30(0x0),
247 fh2JetCoreMethod1C60(0x0),
248 fh2JetCoreMethod2C60(0x0),
249 fh3JetTrackC3060(0x0),
250 fh3JetTrackC20(0x0),
251 fh2AngStructpt1C10(0x0),
252 fh2AngStructpt2C10(0x0),
253 fh2AngStructpt3C10(0x0),
254 fh2AngStructpt4C10(0x0),
255 fh2AngStructpt1C20(0x0),
256 fh2AngStructpt2C20(0x0),
257 fh2AngStructpt3C20(0x0),
258 fh2AngStructpt4C20(0x0),    
259 fh2AngStructpt1C30(0x0),
260 fh2AngStructpt2C30(0x0),
261 fh2AngStructpt3C30(0x0),
262 fh2AngStructpt4C30(0x0),   
263 fh2AngStructpt1C60(0x0),
264 fh2AngStructpt2C60(0x0),
265 fh2AngStructpt3C60(0x0),
266 fh2AngStructpt4C60(0x0),    
267 fh2Ntriggers(0x0),
268 fh2Ntriggers2C10(0x0),
269 fh2Ntriggers2C20(0x0),
270 fh3JetDensity(0x0),
271 fh3JetDensityA4(0x0),
272 fh2RPJetsC10(0x0),
273 fh2RPJetsC20(0x0),
274 fh2RPTC10(0x0),
275 fh2RPTC20(0x0), 
276 fHJetSpec(0x0),
277 fhTTPt(0x0),
278 fHJetPhiCorr(0x0)
279
280  {
281    // Constructor
282
283
284     for(Int_t i=0; i<10; i++) {
285        for(Int_t j=0; j<6; j++) {
286             fTrigBuffer[i][j]=0;
287                 }                               
288     }   
289
290
291
292    fJetBranchName[0] = "";
293    fJetBranchName[1] = "";
294
295    fListJets[0] = new TList;
296    fListJets[1] = new TList;
297
298    DefineOutput(1, TList::Class());
299 }
300
301 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
302 {
303    delete fListJets[0];
304    delete fListJets[1];
305 }
306
307 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
308 {
309    fJetBranchName[0] = branch1;
310    fJetBranchName[1] = branch2;
311 }
312
313 void AliAnalysisTaskJetCore::Init()
314 {
315
316    // check for jet branches
317    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
318       AliError("Jet branch name not set.");
319    }
320
321 }
322
323 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
324 {
325    // Create histograms
326    // Called once
327    OpenFile(1);
328    if(!fOutputList) fOutputList = new TList;
329    fOutputList->SetOwner(kTRUE);
330
331    Bool_t oldStatus = TH1::AddDirectoryStatus();
332    TH1::AddDirectory(kFALSE);
333
334
335    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
336    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
337    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
338    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
339    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
340    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
341    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
342
343      UInt_t entries = 0; // bit coded, see GetDimParams() below 
344      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7; 
345      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
346
347      //change binning in pTtrack
348      Double_t *xPt3 = new Double_t[10];
349      xPt3[0] = 0.;
350      for(int i = 1; i<=9;i++){
351       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
352       else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
353       else xPt3[i] = xPt3[i-1] + 150.; // 18 
354      }
355     fhnDeltaR->SetBinEdges(2,xPt3);
356     delete [] xPt3;
357
358     //change binning in HTI
359      Double_t *xPt4 = new Double_t[14];
360      xPt4[0] = 0.;
361      for(int i = 1; i<=13;i++){
362       if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
363       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
364       else xPt4[i] = xPt4[i-1] + 30.; // 13
365      }
366     fhnDeltaR->SetBinEdges(6,xPt4);
367     delete [] xPt4;
368
369     
370
371
372
373
374    
375      if(fDoEventMixing){    
376      UInt_t cifras = 0; // bit coded, see GetDimParams() below 
377      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<7; 
378      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
379
380     if(fCheckMethods){
381     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
382     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
383     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
384     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
385     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
386     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
387     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
388     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
389
390     fh3JetTrackC3060=new TH3F("JetTrackC3060","",50,0,50,150,0.,150.,35,0.,3.5);
391     fh3JetTrackC20=new TH3F("JetTrackC20","",50,0,50,150,0.,150.,35,0.,3.5);
392     if(fAngStructCloseTracks>0){
393     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
394     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
395     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
396     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
397     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
398     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
399     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
400     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
401     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
402     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
403     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
404     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
405     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
406     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
407     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
408     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
409
410    
411
412     fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
413     fh2Ntriggers2C10=new TH2F("# of triggers2C10","",50,0.,50.,50,0.,50.);
414     fh2Ntriggers2C20=new TH2F("# of triggers2C20","",50,0.,50.,50,0.,50.);
415     fh3JetDensity=new TH3F("Jet density vs mutliplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
416     fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,10,0.,50.);
417     fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
418     fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.); 
419     fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.); 
420     fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);  
421
422
423     
424     
425    fOutputList->Add(fHistEvtSelection);
426
427    fOutputList->Add(fhnDeltaR);
428    
429    fOutputList->Add(fhnMixedEvents);
430          
431      
432   
433       if(fCheckMethods){
434
435       fOutputList->Add(fh2JetCoreMethod1C10);
436       fOutputList->Add(fh2JetCoreMethod2C10);
437       fOutputList->Add(fh2JetCoreMethod1C20);
438       fOutputList->Add(fh2JetCoreMethod2C20);
439       fOutputList->Add(fh2JetCoreMethod1C30);
440       fOutputList->Add(fh2JetCoreMethod2C30);
441       fOutputList->Add(fh2JetCoreMethod1C60);
442       fOutputList->Add(fh2JetCoreMethod2C60);}
443       
444       fOutputList->Add(fh3JetTrackC3060);
445       fOutputList->Add(fh3JetTrackC20);
446      
447             
448      
449
450
451         if(fAngStructCloseTracks>0){
452        fOutputList->Add(fh2AngStructpt1C10);
453        fOutputList->Add(fh2AngStructpt2C10);
454        fOutputList->Add(fh2AngStructpt3C10);
455        fOutputList->Add(fh2AngStructpt4C10); 
456        fOutputList->Add(fh2AngStructpt1C20);
457        fOutputList->Add(fh2AngStructpt2C20);
458        fOutputList->Add(fh2AngStructpt3C20);
459        fOutputList->Add(fh2AngStructpt4C20); 
460        fOutputList->Add(fh2AngStructpt1C30);
461        fOutputList->Add(fh2AngStructpt2C30);
462        fOutputList->Add(fh2AngStructpt3C30);
463        fOutputList->Add(fh2AngStructpt4C30);
464        fOutputList->Add(fh2AngStructpt1C60);
465        fOutputList->Add(fh2AngStructpt2C60);
466        fOutputList->Add(fh2AngStructpt3C60);
467        fOutputList->Add(fh2AngStructpt4C60);}  
468
469
470
471
472  
473         fOutputList->Add(fh2Ntriggers);
474         fOutputList->Add(fh2Ntriggers2C10);
475         fOutputList->Add(fh2Ntriggers2C20); 
476         fOutputList->Add(fh3JetDensity);
477         fOutputList->Add(fh3JetDensityA4);
478         fOutputList->Add(fh2RPJetsC10);
479         fOutputList->Add(fh2RPJetsC20);
480          fOutputList->Add(fh2RPTC10);
481         fOutputList->Add(fh2RPTC20);
482
483         const Int_t dimSpec = 5;
484         const Int_t nBinsSpec[dimSpec]     = {10,100, 140, 50, fNRPBins};
485         const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
486         const Double_t hiBinSpec[dimSpec]  = {100,1, 200, 50, fNRPBins};
487         fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
488         fOutputList->Add(fHJetSpec);  
489
490
491         if(fRunAnaAzimuthalCorrelation)
492           {
493             fhTTPt = new TH2F("fhTTPt","Trigger track p_{T} vs centrality",10,0,100,100,0,100);
494             fOutputList->Add(fhTTPt);
495
496             const Int_t dimCor = 5;
497             const Int_t nBinsCor[dimCor]     = {50, 200, 100,              8,   10};
498             const Double_t lowBinCor[dimCor] = {0,  -50, -0.5*TMath::Pi(), 0,   0};
499             const Double_t hiBinCor[dimCor]  = {50, 150, 1.5*TMath::Pi(),  0.8, 100};
500             fHJetPhiCorr = new THnSparseF("fHJetPhiCorr","TT p_{T} vs jet p_{T} vs dPhi vs area vs centrality",dimCor,nBinsCor,lowBinCor,hiBinCor);
501             fOutputList->Add(fHJetPhiCorr);
502           }
503
504
505         
506      
507    // =========== Switch on Sumw2 for all histos ===========
508    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
509       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
510       if (h1){
511          h1->Sumw2();
512          continue;
513       }
514     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
515       if (hn){
516          hn->Sumw2();
517       }   
518    }
519    TH1::AddDirectory(oldStatus);
520
521    PostData(1, fOutputList);
522 }
523
524 void AliAnalysisTaskJetCore::UserExec(Option_t *)
525 {
526    
527
528    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
529       AliError("Jet branch name not set.");
530       return;
531    }
532
533    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
534    if (!fESD) {
535       AliError("ESD not available");
536       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
537    } 
538       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
539
540        static AliAODEvent* aod = 0;
541        // take all other information from the aod we take the tracks from
542        if(!aod){
543        if(!fESD)aod = fAODIn;
544        else aod = fAODOut;}
545
546    
547  
548     if(fNonStdFile.Length()!=0){
549     // case that we have an AOD extension we need can fetch the jets from the extended output
550     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
551     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
552     if(!fAODExtension){
553       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
554     }}
555     
556
557    // -- event selection --
558    fHistEvtSelection->Fill(1); // number of events before event selection
559
560    // physics selection
561    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
562    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
563          std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
564    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
565       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
566       fHistEvtSelection->Fill(2);
567       PostData(1, fOutputList);
568       return;
569    }
570
571    // vertex selection
572    if(!aod){
573      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
574      fHistEvtSelection->Fill(3);
575       PostData(1, fOutputList);
576    }
577    AliAODVertex* primVtx = aod->GetPrimaryVertex();
578
579    if(!primVtx){
580      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
581      fHistEvtSelection->Fill(3);
582      PostData(1, fOutputList);
583      return;
584    }
585
586    Int_t nTracksPrim = primVtx->GetNContributors();
587    if ((nTracksPrim < fMinContribVtx) ||
588          (primVtx->GetZ() < fVtxZMin) ||
589          (primVtx->GetZ() > fVtxZMax) ){
590       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
591       fHistEvtSelection->Fill(3);
592       PostData(1, fOutputList);
593       return;
594    }
595
596    // event class selection (from jet helper task)
597    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
598    if(fDebug) Printf("Event class %d", eventClass);
599    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
600       fHistEvtSelection->Fill(4);
601       PostData(1, fOutputList);
602       return;
603    }
604
605    // centrality selection
606    AliCentrality *cent = 0x0;
607    Double_t centValue = 0.; 
608    if(fIsPbPb){
609    if(fESD) {cent = fESD->GetCentrality();
610      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
611    else     centValue=aod->GetHeader()->GetCentrality();
612    
613    if(fDebug) printf("centrality: %f\n", centValue);
614       if (centValue < fCentMin || centValue > fCentMax){
615       fHistEvtSelection->Fill(4);
616       PostData(1, fOutputList);
617       return;
618       }}
619
620
621    fHistEvtSelection->Fill(0); 
622    // accepted events  
623    // -- end event selection --
624   
625    // get background
626    AliAODJetEventBackground* externalBackground = 0;
627    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
628       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
629       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
630    }
631    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
632      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
633       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
634    }
635
636     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
637       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
638       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
639     }
640    
641    Float_t rho = 0;
642
643    if(fFlagRandom==0){
644      if(externalBackground)rho = externalBackground->GetBackground(0);}
645    if(fFlagRandom==1){
646       if(externalBackground)rho = externalBackground->GetBackground(2);}
647
648    // fetch jets
649    TClonesArray *aodJets[2];
650    aodJets[0]=0;
651    if(fAODOut&&!aodJets[0]){
652    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); 
653    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
654    if(fAODExtension && !aodJets[0]){ 
655    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
656    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
657      if(fAODIn&&!aodJets[0]){
658    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data())); 
659    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  } 
660
661
662    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
663    //Int_t inord[aodJets[0]->GetEntriesFast()];
664    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
665    //  ptsub[n]=0;
666    //  inord[n]=0;}   
667
668    TList ParticleList;
669    Int_t nT = GetListOfTracks(&ParticleList);
670      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
671       fListJets[iJetType]->Clear();
672       if (!aodJets[iJetType]) continue;
673       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
674       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
675          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
676          if (jet) fListJets[iJetType]->Add(jet);
677          // if(iJetType==0){
678          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
679       }}
680    
681    Double_t etabig=0;
682    Double_t ptbig=0;
683    Double_t areabig=0;
684    Double_t phibig=0.;
685    Double_t etasmall=0;
686    Double_t ptsmall=0;
687    Double_t areasmall=0;
688    Double_t phismall=0.;
689          
690   
691
692    Int_t iCount=0; 
693    Int_t trigJet=-1;
694    Int_t trigBBTrack=-1;
695    Int_t trigInTrack=-1;
696    fRPAngle = aod->GetHeader()->GetEventplane();     
697
698    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);     
699    if(!partback){  
700    PostData(1, fOutputList);
701    return;}
702
703
704    //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
705    //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
706    //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);     
707
708    fh2Ntriggers->Fill(centValue,partback->Pt());
709    Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
710    if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
711    if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
712    Double_t accep=2.*TMath::Pi()*1.8;
713    Int_t injet4=0;
714    Int_t injet=0; 
715    if(fSemigoodCorrect){
716    Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
717    if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){ 
718    PostData(1, fOutputList);
719    return;}
720
721    }
722
723    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
724            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
725            etabig  = jetbig->Eta();
726            phibig  = jetbig->Phi();
727            ptbig   = jetbig->Pt();
728            if(ptbig==0) continue; 
729            Double_t phiBin = RelativePhi(phibig,fRPAngle);       
730            areabig = jetbig->EffectiveAreaCharged();
731            Double_t ptcorr=ptbig-rho*areabig;
732            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
733            if(areabig>=0.07) injet=injet+1;
734            if(areabig>=0.4) injet4=injet4+1;   
735            Double_t dphi=RelativePhi(partback->Phi(),phibig); 
736
737            if(fFlagEtaBkg==1){
738            Double_t etadif= partback->Eta()-etabig;
739            if(TMath::Abs(etadif)<=0.5){             
740           
741            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
742            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
743            if(fFlagEtaBkg==0){
744            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
745            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
746
747
748            if(fFlagJetHadron==0){
749            if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
750            if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
751  
752            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
753
754
755            if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
756            if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
757                    Double_t dismin=100.;
758                    Double_t ptmax=-10.; 
759                    Int_t index1=-1;
760                    Int_t index2=-1;
761           
762                    Float_t phitt=partback->Phi();
763                    if(phitt<0)phitt+=TMath::Pi()*2.; 
764                    Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
765
766                    Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(),phiBintt};
767                   fHJetSpec->Fill(fillspec);
768             
769            
770
771                    if(ptcorr<=0) continue;
772
773                        AliAODTrack* leadtrack=0; 
774                        Int_t ippt=0;
775                        Double_t ppt=-10;
776                        if(fFlagJetHadron==0){   
777                          TRefArray *genTrackList = jetbig->GetRefTracks();
778                          Int_t nTracksGenJet = genTrackList->GetEntriesFast();
779                          AliAODTrack* genTrack;
780                          for(Int_t ir=0; ir<nTracksGenJet; ++ir){
781                            genTrack = (AliAODTrack*)(genTrackList->At(ir));
782                            if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
783                              ippt=ir;}}
784                          leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
785                          if(!leadtrack) continue;
786                        }
787
788
789
790                        AliVParticle* leadtrackb=0;
791                        if(fFlagJetHadron!=0){
792                          Int_t nTb = GetHardestTrackBackToJet(jetbig);
793                          leadtrackb = (AliVParticle*)ParticleList.At(nTb);
794                          if(!leadtrackb) continue;  
795                        }
796
797
798
799
800                        
801                        //store one trigger info                   
802                        if(iCount==0){                        
803                          trigJet=i;
804                          trigBBTrack=nT;
805                          trigInTrack=ippt;
806                          iCount=iCount+1;} 
807                        
808    
809                        if(fCheckMethods){
810                          for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
811                            AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
812                            etasmall  = jetsmall->Eta();
813                            phismall = jetsmall->Phi();
814                            ptsmall   = jetsmall->Pt();
815                            areasmall = jetsmall->EffectiveAreaCharged();
816                            Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
817                            tmpDeltaR=TMath::Sqrt(tmpDeltaR);
818                            //Fraction in the jet core  
819                            if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
820                              index2=j;}  
821                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
822                       index1=j;}} //en of loop over R=0.2 jets
823                 //method1:most concentric jet=core 
824                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
825                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
826                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
827                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
828                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
829                 //method2:hardest contained jet=core   
830                 if(index2!=-1){ 
831                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
832                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
833                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
834                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
835                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
836                   if(centValue<10&&leadtrack) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());                  
837                   if(centValue<20&&leadtrack) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());  
838          if(fDoEventMixing==0 && fFlagOnlyRecoil==0){ 
839          for(int it = 0;it<ParticleList.GetEntries();++it){
840           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
841           Double_t deltaR = jetbig->DeltaR(part);
842           Double_t deltaEta = etabig-part->Eta();
843          
844           Double_t deltaPhi=phibig-part->Phi();
845           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
846           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
847           Double_t pTcont=0;
848           if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
849           if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt(); 
850            Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};  
851            fhnDeltaR->Fill(jetEntries);}
852
853
854           }
855          
856          
857           //end of track loop, we only do it if EM is switched off
858          
859
860
861
862
863        
864
865
866
867    }
868    if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
869    if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
870           //end of jet loop
871
872    //}
873
874
875           if(fDoEventMixing>0){
876             //check before if the trigger exists
877             // fTrigBuffer[i][0] = zvtx
878             // fTrigBuffer[i][1] = phi
879             // fTrigBuffer[i][2] = eta
880             // fTrigBuffer[i][3] = pt_jet
881             // fTrigBuffer[i][4] = pt_trig
882             // fTrigBuffer[i][5]= centrality
883             if(fTindex==10) fTindex=0;
884             if(fTrigBuffer[fTindex][3]>0){
885             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
886             if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){  
887               
888                         for(int it = 0;it<nT;++it){
889                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);         
890                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
891                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
892                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
893                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
894                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
895                         Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};                      
896                         fhnMixedEvents->Fill(triggerEntries);
897                         }
898                         fNevents=fNevents+1;  
899                         if(fNevents==10) fTindex=fTindex+1; 
900             }}}
901
902                if(fTindex==10&&fNevents==10) fCountAgain=0;
903
904                // Copy the triggers from the current event into the buffer.
905                //again, only if the trigger exists:
906                if(fCountAgain==0){
907                 if(trigJet>-1){
908                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                      AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);         
909                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
910                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
911                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
912                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
913                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
914                 fTrigBuffer[fTrigBufferIndex][5] = centValue;
915                 fTrigBufferIndex++;
916                 if(fTrigBufferIndex==9) {fTrigBufferIndex=0; 
917                                          fCountAgain=1;}
918                 }
919                }
920           
921           }
922
923           /////////////////////////////////////////////////////////////////////////////
924           ////////////////////// Rongrong's analysis //////////////////////////////////
925           if(fRunAnaAzimuthalCorrelation)
926             {
927               fhTTPt->Fill(centValue,partback->Pt());
928               for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
929                 {
930                   AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
931                   Double_t jetPt   = jet->Pt();
932                   Double_t jetEta  = jet->Eta();
933                   Double_t jetPhi  = jet->Phi();
934                   if(jetPt==0) continue; 
935                   if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
936                   Double_t jetArea = jet->EffectiveAreaCharged();
937                   Double_t jetPtCorr=jetPt-rho*jetArea;
938                   Double_t dPhi=jetPhi-partback->Phi();
939                   if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
940                   if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
941                   if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
942                   if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
943
944                   Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
945                   fHJetPhiCorr->Fill(fill);
946                 }
947             }
948           /////////////////////////////////////////////////////////////////////////////
949           /////////////////////////////////////////////////////////////////////////////
950
951
952      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
953
954      //tracks up to R=0.8 distant from the jet axis
955    //   if(fAngStructCloseTracks==1){
956    //    TList CloseTrackList;
957    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
958    //    Double_t difR=0.04;
959    //    for(Int_t l=0;l<15;l++){
960    //    Double_t rr=l*0.1+0.1;
961    //     for(int it = 0;it<nn;++it){
962    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
963    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
964    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
965    //         Double_t ptm=part1->Pt();
966    //         Double_t ptn=part2->Pt(); 
967    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
968    //                    Rnm=TMath::Sqrt(Rnm);
969    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
970    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
971    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
972    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
973    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
974    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
975    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
976    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
977    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
978    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
979    //   }
980     
981    //   //only jet constituents
982    //    if(fAngStructCloseTracks==2){
983
984    //    Double_t difR=0.04;
985    //    for(Int_t l=0;l<15;l++){
986    //    Double_t rr=l*0.1+0.1;
987
988     
989    //    AliAODTrack* part1;
990    //    AliAODTrack* part2;
991
992    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
993    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
994           
995              
996
997    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
998    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
999    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
1000    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
1001    //         Double_t ptm=part1->Pt();
1002    //         Double_t ptn=part2->Pt(); 
1003    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1004    //                    Rnm=TMath::Sqrt(Rnm);
1005    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1006    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
1007    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1008    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1009    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1010    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
1011    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1012    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1013    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1014    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1015    // }
1016      // //end loop over R=0.4 jets      
1017      // if(fAngStructCloseTracks>0){
1018      // for(Int_t l=0;l<15;l++){
1019      // Double_t rr=l*0.1+0.1;
1020      //    if(down1[l]!=0){  
1021      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1022      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1023      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1024      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1025      //    if(down2[l]!=0){  
1026      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1027      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1028      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1029      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1030      //    if(down3[l]!=0){  
1031      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1032      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1033      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1034      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1035      //    if(down4[l]!=0){  
1036      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1037      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1038      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1039      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
1040
1041     
1042
1043
1044
1045
1046
1047    PostData(1, fOutputList);
1048 }
1049
1050 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1051 {
1052    // Draw result to the screen
1053    // Called once at the end of the query
1054
1055    if (!GetOutputData(1))
1056    return;
1057 }
1058
1059
1060
1061   
1062
1063
1064 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1065
1066       Int_t iCount = 0;
1067       AliAODEvent *aod = 0;
1068
1069      if(!fESD)aod = fAODIn;
1070      else aod = fAODOut;   
1071
1072      if(!aod)return 0;
1073
1074      Int_t index=-1;
1075      Double_t ptmax=-10;
1076
1077
1078     
1079      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1080       AliAODTrack *tr = aod->GetTrack(it);
1081       Bool_t bGood = false;
1082       if(fFilterType == 0)bGood = true;
1083       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1084       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();    
1085      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1086       if(bGood==false) continue;
1087       if(TMath::Abs(tr->Eta())>0.9)continue;
1088       if(tr->Pt()<0.15)continue;
1089       list->Add(tr);
1090       iCount++;
1091       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1092         if(tr->TestFilterBit(fFilterMaskBestPt)){
1093           if(tr->Pt()>ptmax){ 
1094             ptmax=tr->Pt();     
1095             index=iCount-1;
1096           }
1097         }
1098       }
1099       else{
1100         if(tr->Pt()>ptmax){ 
1101           ptmax=tr->Pt();       
1102           index=iCount-1;
1103         }
1104       }
1105      }
1106   
1107    
1108     // else if (type == kTrackAODMCCharged) {
1109     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1110     // if(!tca)return iCount;
1111     // for(int it = 0;it < tca->GetEntriesFast();++it){
1112     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1113     //   if(!part)continue;
1114     //   if(part->Pt()<0.15)continue;
1115     //   if(!part->IsPhysicalPrimary())continue;
1116     //   if(part->Charge()==0)continue;
1117     //   if(TMath::Abs(part->Eta())>0.9)continue;
1118     //   list->Add(part);
1119     //   iCount++;
1120     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1121     //  index=iCount-1;}}}
1122       return index;
1123  
1124 }
1125
1126    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1127  
1128     AliAODEvent *aod = 0;
1129     if(!fESD)aod = fAODIn;
1130     else aod = fAODOut;     
1131     Int_t index=-1;
1132     Double_t ptmax=-10;
1133     Double_t dphi=0;
1134     Double_t dif=0;
1135     Int_t iCount=0;
1136     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1137       AliAODTrack *tr = aod->GetTrack(it);
1138       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1139       if(TMath::Abs(tr->Eta())>0.9)continue;
1140       if(tr->Pt()<0.15)continue;
1141       iCount=iCount+1;
1142       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
1143       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1144       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1145       index=iCount-1;
1146       dif=dphi;  }}
1147   
1148       return index;
1149
1150    }
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1161
1162     Int_t iCount = 0;
1163       AliAODEvent *aod = 0;
1164      if(!fESD)aod = fAODIn;
1165      else aod = fAODOut;   
1166   
1167       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1168       AliAODTrack *tr = aod->GetTrack(it);
1169       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1170       if(TMath::Abs(tr->Eta())>0.9)continue;
1171       if(tr->Pt()<0.15)continue;
1172       Double_t disR=jetbig->DeltaR(tr);
1173       if(disR>0.8)  continue;
1174       list->Add(tr);
1175       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1176       iCount++;
1177     }
1178   
1179    list->Sort();
1180    return iCount;
1181
1182 }
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1195 {
1196
1197    Int_t nInputTracks = 0;
1198      AliAODEvent *aod = 0;
1199      if(!fESD)aod = fAODIn;
1200      else aod = fAODOut;   
1201    TString jbname(fJetBranchName[1]);
1202    //needs complete event, use jets without background subtraction
1203    for(Int_t i=1; i<=3; ++i){
1204       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1205    }
1206    // use only HI event
1207    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1208    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1209
1210    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1211    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1212    if(!tmpAODjets){
1213       Printf("Jet branch %s not found", jbname.Data());
1214       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1215       return -1;
1216    }
1217    
1218    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1219       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1220       if(!jet) continue;
1221       TRefArray *trackList = jet->GetRefTracks();
1222       Int_t nTracks = trackList->GetEntriesFast();
1223       nInputTracks += nTracks;
1224       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1225    }
1226    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1227
1228    return nInputTracks;  
1229 }
1230
1231
1232
1233 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1234
1235   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1236   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1237   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1238   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1239   double dphi = mphi-vphi;
1240   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1241   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1242   return dphi;//dphi in [-Pi, Pi]
1243 }
1244
1245 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1246 {
1247     Int_t phibin=-1;
1248     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1249     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1250     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1251     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1252     return phibin;
1253 }
1254
1255
1256
1257
1258 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1259 {
1260    // generate new THnSparseF, axes are defined in GetDimParams()
1261
1262    Int_t count = 0;
1263    UInt_t tmp = entries;
1264    while(tmp!=0){
1265       count++;
1266       tmp = tmp &~ -tmp;  // clear lowest bit
1267    }
1268
1269    TString hnTitle(name);
1270    const Int_t dim = count;
1271    Int_t nbins[dim];
1272    Double_t xmin[dim];
1273    Double_t xmax[dim];
1274
1275    Int_t i=0;
1276    Int_t c=0;
1277    while(c<dim && i<32){
1278       if(entries&(1<<i)){
1279       
1280          TString label("");
1281          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1282          hnTitle += Form(";%s",label.Data());
1283          c++;
1284       }
1285       
1286       i++;
1287    }
1288    hnTitle += ";";
1289
1290    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1291 }
1292
1293 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1294 {
1295    // stores label and binning of axis for THnSparse
1296
1297    const Double_t pi = TMath::Pi();
1298    
1299    switch(iEntry){
1300       
1301    case 0:
1302       label = "V0 centrality (%)";
1303      
1304          nbins = 10;
1305          xmin = 0.;
1306          xmax = 100.;
1307          break;
1308       
1309       
1310    case 1:
1311       label = "corrected jet pt";
1312          nbins = 20;
1313          xmin = 0.;
1314          xmax = 200.;
1315           break;
1316       
1317       
1318    case 2:
1319       label = "track pT";
1320      
1321          nbins = 9;
1322          xmin = 0.;
1323          xmax = 150;
1324          break;
1325       
1326       
1327     case 3:
1328       label = "deltaR";
1329       nbins = 15;
1330       xmin = 0.;
1331       xmax = 1.5;
1332       break;
1333
1334
1335
1336    case 4:
1337       label = "deltaEta";
1338       nbins = 8;
1339       xmin = -1.6;
1340       xmax = 1.6;
1341       break;
1342
1343
1344   case 5:
1345       label = "deltaPhi";
1346       nbins = 90;
1347       xmin = -0.5*pi;
1348       xmax = 1.5*pi;
1349       break;   
1350    
1351       
1352         
1353     case 6:
1354       label = "leading track";
1355       nbins = 13;
1356       xmin = 0;
1357       xmax = 50;
1358       break;
1359            
1360      case 7:
1361     
1362       label = "trigger track";
1363       nbins =10;
1364       xmin = 0;
1365       xmax = 50;
1366       break;
1367
1368
1369    
1370   
1371
1372
1373
1374
1375    }
1376
1377 }
1378