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