Coverity
[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    if(!fAOD){
401      PostData(1, fOutputList);
402      if(fDebug) Printf(" No AOD found ... ");
403      return;
404    }
405    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
406    Int_t nTracksPrim = primVtx->GetNContributors();
407    if ((nTracksPrim < fMinContribVtx) ||
408          (primVtx->GetZ() < fVtxZMin) ||
409          (primVtx->GetZ() > fVtxZMax) ){
410       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
411       fHistEvtSelection->Fill(3);
412       PostData(1, fOutputList);
413       return;
414    }
415
416    // event class selection (from jet helper task)
417    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
418    if(fDebug) Printf("Event class %d", eventClass);
419    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
420       fHistEvtSelection->Fill(4);
421       PostData(1, fOutputList);
422       return;
423    }
424
425    // centrality selection
426    AliCentrality *cent = 0x0;
427    Double_t centValue = 0.; 
428    if(fESD) cent = fESD->GetCentrality();
429    if(cent) centValue = cent->GetCentralityPercentile("V0M");
430    if(fDebug) printf("centrality: %f\n", centValue);
431    if (centValue < fCentMin || centValue > fCentMax){
432       fHistEvtSelection->Fill(4);
433       PostData(1, fOutputList);
434       return;
435    }
436
437
438    fHistEvtSelection->Fill(0); 
439    // accepted events  
440    // -- end event selection --
441   
442    // get background
443    AliAODJetEventBackground* externalBackground = 0;
444    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
445       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
446       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
447    }
448    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
449      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
450       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
451    }
452    
453    Float_t rho = 0;
454    if(externalBackground)rho = externalBackground->GetBackground(0);
455
456
457    // fetch jets
458    TClonesArray *aodJets[2];
459    aodJets[0]=0;
460    if(fAOD&&!aodJets[0]){
461    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); 
462    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
463    if(fAODExtension && !aodJets[0]){ 
464      aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
465      aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
466
467    Double_t ptsub[aodJets[0]->GetEntriesFast()];
468    Int_t inord[aodJets[0]->GetEntriesFast()];
469    for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
470      ptsub[n]=0;
471      inord[n]=0;}   
472
473    TList ParticleList;
474    Int_t nT = GetListOfTracks(&ParticleList);
475      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
476       fListJets[iJetType]->Clear();
477       if (!aodJets[iJetType]) continue;
478       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
479       
480    
481       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
482          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
483          if(!jet)continue;
484           fListJets[iJetType]->Add(jet);
485          if(iJetType==0){
486            ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();
487          }}}
488    
489    Double_t etabig=0;
490    Double_t ptbig=0;
491    Double_t areabig=0;
492    Double_t phibig=0.;
493    Double_t etasmall=0;
494    Double_t ptsmall=0;
495    Double_t areasmall=0;
496    //Double_t distr=0.;
497    Double_t phismall=0.;
498    Int_t indexlead=-1;
499    Int_t indexsublead=-1;
500    Int_t indexstop=-1;
501   
502    if(fListJets[0]->GetEntries()>0) TMath::Sort(fListJets[0]->GetEntries(),ptsub,inord);
503
504    for(Int_t jj=0;jj<fListJets[0]->GetEntries();jj++){
505    AliAODJet* jetlead = (AliAODJet*)(fListJets[0]->At(inord[jj]));
506    if(jetlead->Pt()-rho*jetlead->EffectiveAreaCharged()<=0) continue;
507    if((jetlead->Eta()<fJetEtaMin)||(jetlead->Eta()>fJetEtaMax)) continue;
508    indexlead=inord[jj];
509    indexstop=jj;
510    break;}
511    if((indexstop>-1)&&(indexstop+1<fListJets[0]->GetEntries()-1)){
512    for(Int_t k=indexstop+1;k<fListJets[0]->GetEntries();k++){
513    AliAODJet* jetsublead = (AliAODJet*)(fListJets[0]->At(inord[k]));
514    if(jetsublead->Pt()-rho*jetsublead->EffectiveAreaCharged()<=0) continue;
515    if((jetsublead->Eta()<fJetEtaMin)||(jetsublead->Eta()>fJetEtaMax)) continue;
516    indexsublead=inord[k];
517    break;}}
518          
519   
520    Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
521    Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
522    Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
523    Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
524    Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
525    Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
526    Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
527    Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
528     
529   
530
531    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
532            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
533            etabig  = jetbig->Eta();
534            phibig  = jetbig->Phi();
535            ptbig   = jetbig->Pt();
536            if(ptbig==0) continue; 
537            areabig = jetbig->EffectiveAreaCharged();
538            Double_t ptcorr=ptbig-rho*areabig;
539            if(ptcorr<=0) continue;
540            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
541                    Double_t dismin=100.;
542                    Double_t ptmax=-10.; 
543                    Int_t index1=-1;
544                    Int_t index2=-1;
545                    //Double_t fracin=0.;
546                   
547                Int_t point=GetHardestTrackBackToJet(jetbig);    
548                AliVParticle *partback = (AliVParticle*)ParticleList.At(point);                            
549                    if(!partback) continue; 
550                    fh3spectriggered->Fill(centValue,ptcorr,partback->Pt());
551                     Int_t  flaglead=0;
552                        if(i==indexlead) flaglead=1;
553                        if(i==indexsublead) flaglead=2;
554
555                        fh3specleadsublead->Fill(centValue,ptcorr,flaglead);
556
557                        AliAODTrack* leadtrack; 
558                        Int_t ippt=0;
559                        Double_t ppt=-10;   
560
561                         TRefArray *genTrackList = jetbig->GetRefTracks();
562                         Int_t nTracksGenJet = genTrackList->GetEntriesFast();
563                         AliAODTrack* genTrack;
564                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
565                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
566                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
567                          ippt=ir;}}
568                        //Float_t etr=genTrack->Eta();
569                        //Float_t phir=genTrack->Phi();
570                        //distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
571                        //distr=TMath::Sqrt(distr);
572                        //if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
573                         leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
574                         fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
575                        //fhnJetCoreMethod3->Fill(centValue,ptcorr,fracin/ptbig,partback->Pt(),flaglead);
576                        if(fCheckMethods){
577           
578                 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
579                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
580                   etasmall  = jetsmall->Eta();
581                   phismall = jetsmall->Phi();
582                   ptsmall   = jetsmall->Pt();
583                   areasmall = jetsmall->EffectiveAreaCharged();
584                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
585                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
586                      //Fraction in the jet core  
587                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
588                     index2=j;}  
589                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
590                       index1=j;}} //en of loop over R=0.2 jets
591                 //method1:most concentric jet=core 
592                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
593                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
594                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
595                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
596                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
597                 //method2:hardest contained jet=core   
598                 if(index2!=-1){ 
599                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
600                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
601                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
602                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
603                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
604
605                    
606   
607           for(int it = 0;it<nT;++it){
608           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
609           Double_t deltaR = jetbig->DeltaR(part);
610           Double_t deltaEta = etabig-part->Eta();
611           Double_t deltaPhi=phibig-part->Phi();
612           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
613           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
614
615           Double_t jetEntries[9] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,flaglead,leadtrack->Pt(),partback->Pt()};            fhnDeltaR->Fill(jetEntries);
616           }  
617      //end of track loop
618    
619      //fhnSumBkg->Fill(centValue,ptcorr,bkg/jetbig->Pt(),partback->Pt(),flaglead);       
620    
621    
622      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
623
624      //tracks up to R=0.8 distant from the jet axis
625      if(fAngStructCloseTracks==1){
626       TList CloseTrackList;
627       Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
628       Double_t difR=0.04;
629       for(Int_t l=0;l<15;l++){
630       Double_t rr=l*0.1+0.1;
631        for(int it = 0;it<nn;++it){
632            AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
633            for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
634            AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
635            Double_t ptm=part1->Pt();
636            Double_t ptn=part2->Pt();    
637            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
638                       Rnm=TMath::Sqrt(Rnm);
639                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
640                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
641                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
642                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
643                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
644                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
645                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
646                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
647                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
648                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
649      }
650     
651      //only jet constituents
652       if(fAngStructCloseTracks==2){
653
654       Double_t difR=0.04;
655       for(Int_t l=0;l<15;l++){
656       Double_t rr=l*0.1+0.1;
657
658     
659       AliAODTrack* part1;
660       AliAODTrack* part2;
661
662           TRefArray *genTrackListb = jetbig->GetRefTracks();
663           Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
664           
665              
666
667           for(Int_t it=0; it<nTracksGenJetb; ++it){
668              part1 = (AliAODTrack*)(genTrackListb->At(it));
669            for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
670              part2 = (AliAODTrack*)(genTrackListb->At(itu));
671            Double_t ptm=part1->Pt();
672            Double_t ptn=part2->Pt();    
673            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
674                       Rnm=TMath::Sqrt(Rnm);
675                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
676                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
677                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
678                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
679                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
680                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
681                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
682                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
683                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
684                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
685
686
687    
688
689
690
691
692
693
694
695
696
697
698
699
700    }
701
702
703      //end loop over R=0.4 jets 
704      if(fAngStructCloseTracks>0){
705      for(Int_t l=0;l<15;l++){
706      Double_t rr=l*0.1+0.1;
707         if(down1[l]!=0){  
708         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
709         if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
710         if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
711         if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
712         if(down2[l]!=0){  
713         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
714         if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
715         if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
716         if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
717         if(down3[l]!=0){  
718         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
719         if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
720         if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
721         if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
722         if(down4[l]!=0){  
723         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
724         if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
725         if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
726         if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
727
728     
729
730
731
732
733
734    PostData(1, fOutputList);
735 }
736
737 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
738 {
739    // Draw result to the screen
740    // Called once at the end of the query
741
742    if (!GetOutputData(1))
743    return;
744 }
745
746
747
748
749
750
751
752
753
754
755
756 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
757
758     Int_t iCount = 0;
759  
760     
761     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
762       AliAODTrack *tr = fAOD->GetTrack(it);
763       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
764       if(TMath::Abs(tr->Eta())>0.9)continue;
765       if(tr->Pt()<0.15)continue;
766       list->Add(tr);
767       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
768       iCount++;
769     }
770   
771    list->Sort();
772   return iCount;
773  
774 }
775
776    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
777
778    
779     Int_t index=-1;
780     Double_t ptmax=-10;
781     Double_t dphi=0;
782     Double_t dif=0;
783     Int_t iCount=0;
784     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
785       AliAODTrack *tr = fAOD->GetTrack(it);
786       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
787       if(TMath::Abs(tr->Eta())>0.9)continue;
788       if(tr->Pt()<0.15)continue;
789       iCount=iCount+1;
790       dphi=RelativePhi(tr->Phi(),jetbig->Phi());   
791       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
792       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
793         index=iCount;
794         dif=dphi;  }}
795   
796       return index;
797
798    }
799
800
801
802
803
804
805
806
807
808  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
809
810     Int_t iCount = 0;
811  
812   
813     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
814       AliAODTrack *tr = fAOD->GetTrack(it);
815       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
816       if(TMath::Abs(tr->Eta())>0.9)continue;
817       if(tr->Pt()<0.15)continue;
818       Double_t disR=jetbig->DeltaR(tr);
819       if(disR>0.8)  continue;
820       list->Add(tr);
821       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
822       iCount++;
823     }
824   
825    list->Sort();
826    return iCount;
827
828 }
829
830
831
832
833
834
835
836
837
838
839
840 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
841 {
842
843    Int_t nInputTracks = 0;
844
845    TString jbname(fJetBranchName[1]);
846    //needs complete event, use jets without background subtraction
847    for(Int_t i=1; i<=3; ++i){
848       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
849    }
850    // use only HI event
851    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
852    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
853
854    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
855    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
856    if(!tmpAODjets){
857       Printf("Jet branch %s not found", jbname.Data());
858       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
859       return -1;
860    }
861    
862    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
863       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
864       if(!jet) continue;
865       TRefArray *trackList = jet->GetRefTracks();
866       Int_t nTracks = trackList->GetEntriesFast();
867       nInputTracks += nTracks;
868       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
869    }
870    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
871
872    return nInputTracks;  
873 }
874
875
876
877 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
878
879   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
880   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
881   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
882   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
883   double dphi = mphi-vphi;
884   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
885   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
886   return dphi;//dphi in [-Pi, Pi]
887 }
888
889
890
891 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
892 {
893    // generate new THnSparseF, axes are defined in GetDimParams()
894
895    Int_t count = 0;
896    UInt_t tmp = entries;
897    while(tmp!=0){
898       count++;
899       tmp = tmp &~ -tmp;  // clear lowest bit
900    }
901
902    TString hnTitle(name);
903    const Int_t dim = count;
904    Int_t nbins[dim];
905    Double_t xmin[dim];
906    Double_t xmax[dim];
907
908    Int_t i=0;
909    Int_t c=0;
910    while(c<dim && i<32){
911       if(entries&(1<<i)){
912       
913          TString label("");
914          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
915          hnTitle += Form(";%s",label.Data());
916          c++;
917       }
918       
919       i++;
920    }
921    hnTitle += ";";
922
923    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
924 }
925
926 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
927 {
928    // stores label and binning of axis for THnSparse
929
930    const Double_t pi = TMath::Pi();
931    
932    switch(iEntry){
933       
934    case 0:
935       label = "V0 centrality (%)";
936      
937          nbins = 10;
938          xmin = 0.;
939          xmax = 100.;
940          break;
941       
942       
943    case 1:
944       label = "corrected jet pt";
945          nbins = 50;
946          xmin = 0.;
947          xmax = 200.;
948           break;
949       
950       
951    case 2:
952       label = "track pT";
953      
954       nbins = 50;
955          xmin = 0.;
956          xmax = 50;
957          break;
958       
959       
960    case 3:
961       label = "deltaR";
962       nbins = 15;
963       xmin = 0.;
964       xmax = 1.5;
965       break;
966       
967    case 4:
968       label = "deltaEta";
969       nbins = 30;
970       xmin = -1.5;
971       xmax = 1.5;
972       break;
973
974
975   case 5:
976       label = "deltaPhi";
977       nbins = 90;
978       xmin = -0.5*pi;
979       xmax = 1.5*pi;
980       break;   
981    
982       
983    case 6:
984       label="flagleadname";
985       nbins=3;
986       xmin=0;
987       xmax=3;
988       break;
989
990     case 7:
991     
992       label = "leading track";
993       nbins = 50;
994       xmin = 0;
995       xmax = 50;
996       break;
997            
998      case 8:
999     
1000       label = "trigger track";
1001       nbins =50;
1002       xmin = 0;
1003       xmax = 50;
1004       break;
1005       
1006         
1007    }
1008
1009 }
1010