94ca97e93f6e16305d1fe772bb3e6b92fbfe7f3d
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
1
2 // ******************************************
3 // This task computes several jet observables like 
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and 
6 // correlation strength distribution of particles inside jets.    
7 // Author: lcunquei@cern.ch
8 // *******************************************
9
10
11 /**************************************************************************
12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13  *                                                                        *
14  * Author: The ALICE Off-line Project.                                    *
15  * Contributors are mentioned in the code where appropriate.              *
16  *                                                                        *
17  * Permission to use, copy, modify and distribute this software and its   *
18  * documentation strictly for non-commercial purposes is hereby granted   *
19  * without fee, provided that the above copyright notice appears in all   *
20  * copies and that both the copyright notice and this permission notice   *
21  * appear in the supporting documentation. The authors make no claims     *
22  * about the suitability of this software for any purpose. It is          *
23  * provided "as is" without express or implied warranty.                  *
24  **************************************************************************/
25
26
27 #include "TChain.h"
28 #include "TTree.h"
29 #include "TMath.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH3F.h"
33 #include "THnSparse.h"
34 #include "TCanvas.h"
35
36 #include "AliLog.h"
37
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliVEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliCentrality.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliInputEventHandler.h"
47 #include "AliAODJetEventBackground.h"
48 #include "AliAnalysisTaskFastEmbedding.h"
49 #include "AliAODEvent.h"
50 #include "AliAODHandler.h"
51 #include "AliAODJet.h"
52
53 #include "AliAnalysisTaskJetCore.h"
54
55 ClassImp(AliAnalysisTaskJetCore)
56
57 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
58 AliAnalysisTaskSE(),
59 fESD(0x0),
60 fAOD(0x0),
61 fAODExtension(0x0),
62 fBackgroundBranch(""),
63 fNonStdFile(""),
64 fIsPbPb(kTRUE),
65 fOfflineTrgMask(AliVEvent::kAny),
66 fMinContribVtx(1),
67 fVtxZMin(-10.),
68 fVtxZMax(10.),
69 fEvtClassMin(0),
70 fEvtClassMax(4),
71 fFilterMask(0),
72 fRadioFrac(0.2),
73 fMinDist(0.1),
74 fCentMin(0.),
75 fCentMax(100.),
76 fNInputTracksMin(0),
77 fNInputTracksMax(-1),
78 fAngStructCloseTracks(0),
79 fCheckMethods(0),
80 fJetEtaMin(-.5),
81 fJetEtaMax(.5),
82 fJetPtMin(20.),
83 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
84 fJetPtFractionMin(0.5),
85 fNMatchJets(4),
86 fMatchMaxDist(0.8),
87 fKeepJets(kFALSE),
88 fkNbranches(2),
89 fkEvtClasses(12),
90 fOutputList(0x0),
91 fbEvent(kTRUE),
92 fHistEvtSelection(0x0),
93 fhnDeltaR(0x0),
94 fhnSumBkg(0x0),  
95 fhnJetCoreMethod3(0x0),
96 fh2JetCoreMethod1C10(0x0),
97 fh2JetCoreMethod2C10(0x0),
98 fh2JetCoreMethod1C20(0x0),
99 fh2JetCoreMethod2C20(0x0),
100 fh2JetCoreMethod1C30(0x0),
101 fh2JetCoreMethod2C30(0x0),
102 fh2JetCoreMethod1C60(0x0),
103 fh2JetCoreMethod2C60(0x0),
104 fh2AngStructpt1C10(0x0),
105 fh2AngStructpt2C10(0x0),
106 fh2AngStructpt3C10(0x0),
107 fh2AngStructpt4C10(0x0),
108 fh2AngStructpt1C20(0x0),
109 fh2AngStructpt2C20(0x0),
110 fh2AngStructpt3C20(0x0),
111 fh2AngStructpt4C20(0x0),    
112 fh2AngStructpt1C30(0x0),
113 fh2AngStructpt2C30(0x0),
114 fh2AngStructpt3C30(0x0),
115 fh2AngStructpt4C30(0x0),   
116 fh2AngStructpt1C60(0x0),
117 fh2AngStructpt2C60(0x0),
118 fh2AngStructpt3C60(0x0),
119 fh2AngStructpt4C60(0x0),
120 fh3spectriggered(0x0),
121 fh3specbiased(0x0),
122 fh3specleadsublead(0x0)
123  
124 {
125    // default Constructor
126
127    fJetBranchName[0] = "";
128    fJetBranchName[1] = "";
129
130    fListJets[0] = new TList;
131    fListJets[1] = new TList;
132 }
133
134 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
135 AliAnalysisTaskSE(name),
136 fESD(0x0),
137 fAOD(0x0),
138 fAODExtension(0x0),
139 fBackgroundBranch(""),
140 fNonStdFile(""),
141 fIsPbPb(kTRUE),
142 fOfflineTrgMask(AliVEvent::kAny),
143 fMinContribVtx(1),
144 fVtxZMin(-10.),
145 fVtxZMax(10.),
146 fEvtClassMin(0),
147 fEvtClassMax(4),
148 fFilterMask(0),
149 fRadioFrac(0.2),
150 fMinDist(0.1),
151 fCentMin(0.),
152 fCentMax(100.),
153 fNInputTracksMin(0),
154 fNInputTracksMax(-1),
155 fAngStructCloseTracks(0),
156 fCheckMethods(0),
157 fJetEtaMin(-.5),
158 fJetEtaMax(.5),
159 fJetPtMin(20.),
160 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
161 fJetPtFractionMin(0.5),
162 fNMatchJets(4),
163 fMatchMaxDist(0.8),
164 fKeepJets(kFALSE),
165 fkNbranches(2),
166 fkEvtClasses(12),
167 fOutputList(0x0),
168 fbEvent(kTRUE),
169 fHistEvtSelection(0x0),
170 fhnDeltaR(0x0),
171 fhnSumBkg(0x0),  
172 fhnJetCoreMethod3(0x0),
173 fh2JetCoreMethod1C10(0x0),
174 fh2JetCoreMethod2C10(0x0),
175 fh2JetCoreMethod1C20(0x0),
176 fh2JetCoreMethod2C20(0x0),
177 fh2JetCoreMethod1C30(0x0),
178 fh2JetCoreMethod2C30(0x0),
179 fh2JetCoreMethod1C60(0x0),
180 fh2JetCoreMethod2C60(0x0),
181 fh2AngStructpt1C10(0x0),
182 fh2AngStructpt2C10(0x0),
183 fh2AngStructpt3C10(0x0),
184 fh2AngStructpt4C10(0x0),
185 fh2AngStructpt1C20(0x0),
186 fh2AngStructpt2C20(0x0),
187 fh2AngStructpt3C20(0x0),
188 fh2AngStructpt4C20(0x0),    
189 fh2AngStructpt1C30(0x0),
190 fh2AngStructpt2C30(0x0),
191 fh2AngStructpt3C30(0x0),
192 fh2AngStructpt4C30(0x0),   
193 fh2AngStructpt1C60(0x0),
194 fh2AngStructpt2C60(0x0),
195 fh2AngStructpt3C60(0x0),
196 fh2AngStructpt4C60(0x0),    
197 fh3spectriggered(0x0),
198 fh3specbiased(0x0),
199 fh3specleadsublead(0x0)
200  {
201    // Constructor
202
203    fJetBranchName[0] = "";
204    fJetBranchName[1] = "";
205
206    fListJets[0] = new TList;
207    fListJets[1] = new TList;
208
209    DefineOutput(1, TList::Class());
210 }
211
212 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
213 {
214    delete fListJets[0];
215    delete fListJets[1];
216 }
217
218 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
219 {
220    fJetBranchName[0] = branch1;
221    fJetBranchName[1] = branch2;
222 }
223
224 void AliAnalysisTaskJetCore::Init()
225 {
226
227    // check for jet branches
228    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
229       AliError("Jet branch name not set.");
230    }
231
232 }
233
234 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
235 {
236    // Create histograms
237    // Called once
238    OpenFile(1);
239    if(!fOutputList) fOutputList = new TList;
240    fOutputList->SetOwner(kTRUE);
241
242    Bool_t oldStatus = TH1::AddDirectoryStatus();
243    TH1::AddDirectory(kFALSE);
244
245
246    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
247    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
248    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
249    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
250    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
251    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
252    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
253
254     UInt_t entries = 0; // bit coded, see GetDimParams() below 
255     
256      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7|1<<8; 
257      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
258      // fhnSumBkg = NewTHnSparseF("fhnDeltaR", entries);
259
260     if(fCheckMethods){
261       entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5  ; 
262      fhnJetCoreMethod3 = NewTHnSparseF("fhnEvent", entries);
263
264     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
265     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
266     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
267     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
268     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
269     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
270     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
271     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
272
273    
274     if(fAngStructCloseTracks>0){
275     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
276     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
277     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
278     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
279     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
280     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
281     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
282     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
283     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
284     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
285     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
286     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
287     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
288     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
289     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
290     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
291     fh3spectriggered = new TH3F("Triggered spectrum","",10,0,100,50,0.,200,50,0.,50.);
292     fh3specbiased = new TH3F("Biased spectrum","",10,0,100,50,0.,200.,50,0.,50.);
293     fh3specleadsublead = new TH3F("Leading/subleading spectrum","",10,0,100,50,0.,200.,3,0,3);
294
295    fOutputList->Add(fHistEvtSelection);
296
297    fOutputList->Add(fhnDeltaR);
298    
299    //fOutputList->Add(fhnSumBkg);
300          
301      
302   
303       if(fCheckMethods){
304       fOutputList->Add(fhnJetCoreMethod3);
305       fOutputList->Add(fh2JetCoreMethod1C10);
306       fOutputList->Add(fh2JetCoreMethod2C10);
307       fOutputList->Add(fh2JetCoreMethod1C20);
308       fOutputList->Add(fh2JetCoreMethod2C20);
309       fOutputList->Add(fh2JetCoreMethod1C30);
310       fOutputList->Add(fh2JetCoreMethod2C30);
311       fOutputList->Add(fh2JetCoreMethod1C60);
312       fOutputList->Add(fh2JetCoreMethod2C60);}
313       
314       
315      
316
317
318         if(fAngStructCloseTracks>0){
319        fOutputList->Add(fh2AngStructpt1C10);
320        fOutputList->Add(fh2AngStructpt2C10);
321        fOutputList->Add(fh2AngStructpt3C10);
322        fOutputList->Add(fh2AngStructpt4C10); 
323        fOutputList->Add(fh2AngStructpt1C20);
324        fOutputList->Add(fh2AngStructpt2C20);
325        fOutputList->Add(fh2AngStructpt3C20);
326        fOutputList->Add(fh2AngStructpt4C20); 
327        fOutputList->Add(fh2AngStructpt1C30);
328        fOutputList->Add(fh2AngStructpt2C30);
329        fOutputList->Add(fh2AngStructpt3C30);
330        fOutputList->Add(fh2AngStructpt4C30);
331        fOutputList->Add(fh2AngStructpt1C60);
332        fOutputList->Add(fh2AngStructpt2C60);
333        fOutputList->Add(fh2AngStructpt3C60);
334        fOutputList->Add(fh2AngStructpt4C60);}  
335        fOutputList->Add(fh3spectriggered);
336        fOutputList->Add(fh3specbiased);
337        fOutputList->Add(fh3specleadsublead);
338
339    // =========== Switch on Sumw2 for all histos ===========
340    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
341       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
342       if (h1){
343          h1->Sumw2();
344          continue;
345       }
346     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
347       if (hn){
348          hn->Sumw2();
349       }   
350    }
351    TH1::AddDirectory(oldStatus);
352
353    PostData(1, fOutputList);
354 }
355
356 void AliAnalysisTaskJetCore::UserExec(Option_t *)
357 {
358    
359
360    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
361       AliError("Jet branch name not set.");
362       return;
363    }
364
365    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
366    if (!fESD) {
367       AliError("ESD not available");
368       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
369    } else {
370       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
371    }
372  
373     if(fNonStdFile.Length()!=0){
374     // case that we have an AOD extension we need can fetch the jets from the extended output
375     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
376     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
377     if(!fAODExtension){
378       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
379     }}
380     
381
382
383
384
385    // -- event selection --
386    fHistEvtSelection->Fill(1); // number of events before event selection
387
388    // physics selection
389    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
390    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
391    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
392    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
393       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
394       fHistEvtSelection->Fill(2);
395       PostData(1, fOutputList);
396       return;
397    }
398
399    // vertex selection
400    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
401    Int_t nTracksPrim = primVtx->GetNContributors();
402    if ((nTracksPrim < fMinContribVtx) ||
403          (primVtx->GetZ() < fVtxZMin) ||
404          (primVtx->GetZ() > fVtxZMax) ){
405       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
406       fHistEvtSelection->Fill(3);
407       PostData(1, fOutputList);
408       return;
409    }
410
411    // event class selection (from jet helper task)
412    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
413    if(fDebug) Printf("Event class %d", eventClass);
414    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
415       fHistEvtSelection->Fill(4);
416       PostData(1, fOutputList);
417       return;
418    }
419
420    // centrality selection
421    AliCentrality *cent = 0x0;
422    Double_t centValue = 0.; 
423    if(fESD) cent = fESD->GetCentrality();
424    if(cent) centValue = cent->GetCentralityPercentile("V0M");
425    if(fDebug) printf("centrality: %f\n", centValue);
426    if (centValue < fCentMin || centValue > fCentMax){
427       fHistEvtSelection->Fill(4);
428       PostData(1, fOutputList);
429       return;
430    }
431
432
433    fHistEvtSelection->Fill(0); 
434    // accepted events  
435    // -- end event selection --
436   
437    // get background
438    AliAODJetEventBackground* externalBackground = 0;
439    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
440       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
441       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
442    }
443    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
444      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
445       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
446    }
447    
448    Float_t rho = 0;
449    if(externalBackground)rho = externalBackground->GetBackground(0);
450
451
452    // fetch jets
453    TClonesArray *aodJets[2];
454    aodJets[0]=0;
455    if(fAOD&&!aodJets[0]){
456    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); 
457    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
458    if(fAODExtension && !aodJets[0]){ 
459      aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
460      aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
461
462    Double_t ptsub[aodJets[0]->GetEntriesFast()];
463    Int_t inord[aodJets[0]->GetEntriesFast()];
464    for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
465      ptsub[n]=0;
466      inord[n]=0;}   
467
468    TList ParticleList;
469    Int_t nT = GetListOfTracks(&ParticleList);
470      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
471       fListJets[iJetType]->Clear();
472       if (!aodJets[iJetType]) continue;
473
474       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
475       
476    
477       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
478          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
479          if (jet) fListJets[iJetType]->Add(jet);
480          if(iJetType==0){
481            ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();
482          }}}
483    
484    Double_t etabig=0;
485    Double_t ptbig=0;
486    Double_t areabig=0;
487    Double_t phibig=0.;
488    Double_t etasmall=0;
489    Double_t ptsmall=0;
490    Double_t areasmall=0;
491    //Double_t distr=0.;
492    Double_t phismall=0.;
493    Int_t indexlead=-1;
494    Int_t indexsublead=-1;
495    Int_t indexstop=-1;
496   
497    if(fListJets[0]->GetEntries()>0) TMath::Sort(fListJets[0]->GetEntries(),ptsub,inord);
498
499    for(Int_t jj=0;jj<fListJets[0]->GetEntries();jj++){
500    AliAODJet* jetlead = (AliAODJet*)(fListJets[0]->At(inord[jj]));
501    if(jetlead->Pt()-rho*jetlead->EffectiveAreaCharged()<=0) continue;
502    if((jetlead->Eta()<fJetEtaMin)||(jetlead->Eta()>fJetEtaMax)) continue;
503    indexlead=inord[jj];
504    indexstop=jj;
505    break;}
506    if((indexstop>-1)&&(indexstop+1<fListJets[0]->GetEntries()-1)){
507    for(Int_t k=indexstop+1;k<fListJets[0]->GetEntries();k++){
508    AliAODJet* jetsublead = (AliAODJet*)(fListJets[0]->At(inord[k]));
509    if(jetsublead->Pt()-rho*jetsublead->EffectiveAreaCharged()<=0) continue;
510    if((jetsublead->Eta()<fJetEtaMin)||(jetsublead->Eta()>fJetEtaMax)) continue;
511    indexsublead=inord[k];
512    break;}}
513          
514   
515    Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
516    Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
517    Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
518    Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
519    Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
520    Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
521    Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
522    Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
523     
524   
525
526    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
527            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
528            etabig  = jetbig->Eta();
529            phibig  = jetbig->Phi();
530            ptbig   = jetbig->Pt();
531            if(ptbig==0) continue; 
532            areabig = jetbig->EffectiveAreaCharged();
533            Double_t ptcorr=ptbig-rho*areabig;
534            if(ptcorr<=0) continue;
535            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
536                    Double_t dismin=100.;
537                    Double_t ptmax=-10.; 
538                    Int_t index1=-1;
539                    Int_t index2=-1;
540                    //Double_t fracin=0.;
541                   
542                Int_t point=GetHardestTrackBackToJet(jetbig);    
543                AliVParticle *partback = (AliVParticle*)ParticleList.At(point);                            
544                    if(!partback) continue; 
545                    fh3spectriggered->Fill(centValue,ptcorr,partback->Pt());
546                     Int_t  flaglead=0;
547                        if(i==indexlead) flaglead=1;
548                        if(i==indexsublead) flaglead=2;
549
550                        fh3specleadsublead->Fill(centValue,ptcorr,flaglead);
551
552                        AliAODTrack* leadtrack; 
553                        Int_t ippt=0;
554                        Double_t ppt=-10;   
555
556                         TRefArray *genTrackList = jetbig->GetRefTracks();
557                         Int_t nTracksGenJet = genTrackList->GetEntriesFast();
558                         AliAODTrack* genTrack;
559                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
560                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
561                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
562                          ippt=ir;}}
563                        //Float_t etr=genTrack->Eta();
564                        //Float_t phir=genTrack->Phi();
565                        //distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
566                        //distr=TMath::Sqrt(distr);
567                        //if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
568                         leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
569                         fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
570                        //fhnJetCoreMethod3->Fill(centValue,ptcorr,fracin/ptbig,partback->Pt(),flaglead);
571                        if(fCheckMethods){
572           
573                 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
574                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
575                   etasmall  = jetsmall->Eta();
576                   phismall = jetsmall->Phi();
577                   ptsmall   = jetsmall->Pt();
578                   areasmall = jetsmall->EffectiveAreaCharged();
579                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
580                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
581                      //Fraction in the jet core  
582                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
583                     index2=j;}  
584                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
585                       index1=j;}} //en of loop over R=0.2 jets
586                 //method1:most concentric jet=core 
587                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
588                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
589                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
590                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
591                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
592                 //method2:hardest contained jet=core   
593                 if(index2!=-1){ 
594                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
595                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
596                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
597                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
598                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
599
600                    
601   
602           for(int it = 0;it<nT;++it){
603           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
604           Double_t deltaR = jetbig->DeltaR(part);
605           Double_t deltaEta = etabig-part->Eta();
606           Double_t deltaPhi=phibig-part->Phi();
607           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
608           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
609
610           Double_t jetEntries[9] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,flaglead,leadtrack->Pt(),partback->Pt()};            fhnDeltaR->Fill(jetEntries);
611           }  
612      //end of track loop
613    
614      //fhnSumBkg->Fill(centValue,ptcorr,bkg/jetbig->Pt(),partback->Pt(),flaglead);       
615    
616    
617      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
618
619      //tracks up to R=0.8 distant from the jet axis
620      if(fAngStructCloseTracks==1){
621       TList CloseTrackList;
622       Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
623       Double_t difR=0.04;
624       for(Int_t l=0;l<15;l++){
625       Double_t rr=l*0.1+0.1;
626        for(int it = 0;it<nn;++it){
627            AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
628            for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
629            AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
630            Double_t ptm=part1->Pt();
631            Double_t ptn=part2->Pt();    
632            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
633                       Rnm=TMath::Sqrt(Rnm);
634                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
635                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
636                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
637                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
638                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
639                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
640                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
641                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
642                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
643                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
644      }
645     
646      //only jet constituents
647       if(fAngStructCloseTracks==2){
648
649       Double_t difR=0.04;
650       for(Int_t l=0;l<15;l++){
651       Double_t rr=l*0.1+0.1;
652
653     
654       AliAODTrack* part1;
655       AliAODTrack* part2;
656
657           TRefArray *genTrackListb = jetbig->GetRefTracks();
658           Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
659           
660              
661
662           for(Int_t it=0; it<nTracksGenJetb; ++it){
663              part1 = (AliAODTrack*)(genTrackListb->At(it));
664            for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
665              part2 = (AliAODTrack*)(genTrackListb->At(itu));
666            Double_t ptm=part1->Pt();
667            Double_t ptn=part2->Pt();    
668            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
669                       Rnm=TMath::Sqrt(Rnm);
670                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
671                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
672                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
673                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
674                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
675                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
676                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
677                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
678                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
679                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
680
681
682    
683
684
685
686
687
688
689
690
691
692
693
694
695    }
696
697
698      //end loop over R=0.4 jets 
699      if(fAngStructCloseTracks>0){
700      for(Int_t l=0;l<15;l++){
701      Double_t rr=l*0.1+0.1;
702         if(down1[l]!=0){  
703         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
704         if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
705         if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
706         if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
707         if(down2[l]!=0){  
708         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
709         if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
710         if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
711         if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
712         if(down3[l]!=0){  
713         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
714         if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
715         if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
716         if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
717         if(down4[l]!=0){  
718         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
719         if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
720         if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
721         if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
722
723     
724
725
726
727
728
729    PostData(1, fOutputList);
730 }
731
732 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
733 {
734    // Draw result to the screen
735    // Called once at the end of the query
736
737    if (!GetOutputData(1))
738    return;
739 }
740
741
742
743
744
745
746
747
748
749
750
751 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
752
753     Int_t iCount = 0;
754  
755     
756     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
757       AliAODTrack *tr = fAOD->GetTrack(it);
758       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
759       if(TMath::Abs(tr->Eta())>0.9)continue;
760       if(tr->Pt()<0.15)continue;
761       list->Add(tr);
762       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
763       iCount++;
764     }
765   
766    list->Sort();
767   return iCount;
768  
769 }
770
771    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
772
773    
774     Int_t index=-1;
775     Double_t ptmax=-10;
776     Double_t dphi=0;
777     Double_t dif=0;
778     Int_t iCount=0;
779     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
780       AliAODTrack *tr = fAOD->GetTrack(it);
781       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
782       if(TMath::Abs(tr->Eta())>0.9)continue;
783       if(tr->Pt()<0.15)continue;
784       iCount=iCount+1;
785       dphi=RelativePhi(tr->Phi(),jetbig->Phi());   
786       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
787       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
788         index=iCount;
789         dif=dphi;  }}
790   
791       return index;
792
793    }
794
795
796
797
798
799
800
801
802
803  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
804
805     Int_t iCount = 0;
806  
807   
808     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
809       AliAODTrack *tr = fAOD->GetTrack(it);
810       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
811       if(TMath::Abs(tr->Eta())>0.9)continue;
812       if(tr->Pt()<0.15)continue;
813       Double_t disR=jetbig->DeltaR(tr);
814       if(disR>0.8)  continue;
815       list->Add(tr);
816       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
817       iCount++;
818     }
819   
820    list->Sort();
821    return iCount;
822
823 }
824
825
826
827
828
829
830
831
832
833
834
835 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
836 {
837
838    Int_t nInputTracks = 0;
839
840    TString jbname(fJetBranchName[1]);
841    //needs complete event, use jets without background subtraction
842    for(Int_t i=1; i<=3; ++i){
843       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
844    }
845    // use only HI event
846    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
847    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
848
849    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
850    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
851    if(!tmpAODjets){
852       Printf("Jet branch %s not found", jbname.Data());
853       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
854       return -1;
855    }
856    
857    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
858       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
859       if(!jet) continue;
860       TRefArray *trackList = jet->GetRefTracks();
861       Int_t nTracks = trackList->GetEntriesFast();
862       nInputTracks += nTracks;
863       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
864    }
865    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
866
867    return nInputTracks;  
868 }
869
870
871
872 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
873
874   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
875   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
876   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
877   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
878   double dphi = mphi-vphi;
879   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
880   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
881   return dphi;//dphi in [-Pi, Pi]
882 }
883
884
885
886 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
887 {
888    // generate new THnSparseF, axes are defined in GetDimParams()
889
890    Int_t count = 0;
891    UInt_t tmp = entries;
892    while(tmp!=0){
893       count++;
894       tmp = tmp &~ -tmp;  // clear lowest bit
895    }
896
897    TString hnTitle(name);
898    const Int_t dim = count;
899    Int_t nbins[dim];
900    Double_t xmin[dim];
901    Double_t xmax[dim];
902
903    Int_t i=0;
904    Int_t c=0;
905    while(c<dim && i<32){
906       if(entries&(1<<i)){
907       
908          TString label("");
909          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
910          hnTitle += Form(";%s",label.Data());
911          c++;
912       }
913       
914       i++;
915    }
916    hnTitle += ";";
917
918    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
919 }
920
921 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
922 {
923    // stores label and binning of axis for THnSparse
924
925    const Double_t pi = TMath::Pi();
926    
927    switch(iEntry){
928       
929    case 0:
930       label = "V0 centrality (%)";
931      
932          nbins = 10;
933          xmin = 0.;
934          xmax = 100.;
935          break;
936       
937       
938    case 1:
939       label = "corrected jet pt";
940          nbins = 50;
941          xmin = 0.;
942          xmax = 200.;
943           break;
944       
945       
946    case 2:
947       label = "track pT";
948      
949       nbins = 50;
950          xmin = 0.;
951          xmax = 50;
952          break;
953       
954       
955    case 3:
956       label = "deltaR";
957       nbins = 15;
958       xmin = 0.;
959       xmax = 1.5;
960       break;
961       
962    case 4:
963       label = "deltaEta";
964       nbins = 30;
965       xmin = -1.5;
966       xmax = 1.5;
967       break;
968
969
970   case 5:
971       label = "deltaPhi";
972       nbins = 90;
973       xmin = -0.5*pi;
974       xmax = 1.5*pi;
975       break;   
976    
977       
978    case 6:
979       label="flagleadname";
980       nbins=3;
981       xmin=0;
982       xmax=3;
983
984
985
986       
987     case 7:
988     
989       label = "leading track";
990       nbins = 50;
991       xmin = 0;
992       xmax = 50;
993       break;
994            
995      case 8:
996     
997       label = "trigger track";
998       nbins =50;
999       xmin = 0;
1000       xmax = 50;
1001       break;
1002       
1003         
1004    }
1005
1006 }
1007