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