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