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