added jet trigger selection and jet hadron correlation
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
CommitLineData
a9e585a7 1
568f8fa2 2// ******************************************
3// This task computes several jet observables like
4// the fraction of energy in inner and outer coronnas,
da93bb11 5// jet-track correlations,triggered jet shapes and
568f8fa2 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
75bf77e3 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"
75bf77e3 49#include "AliAODEvent.h"
ea693273 50#include "AliAODHandler.h"
75bf77e3 51#include "AliAODJet.h"
52
53#include "AliAnalysisTaskJetCore.h"
54
55ClassImp(AliAnalysisTaskJetCore)
56
57AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
58AliAnalysisTaskSE(),
59fESD(0x0),
60fAOD(0x0),
ea693273 61fAODExtension(0x0),
75bf77e3 62fBackgroundBranch(""),
ea693273 63fNonStdFile(""),
75bf77e3 64fIsPbPb(kTRUE),
65fOfflineTrgMask(AliVEvent::kAny),
66fMinContribVtx(1),
a9e585a7 67fVtxZMin(-10.),
68fVtxZMax(10.),
75bf77e3 69fEvtClassMin(0),
70fEvtClassMax(4),
8b47ec90 71fFilterMask(0),
ea693273 72fRadioFrac(0.2),
73fMinDist(0.1),
75bf77e3 74fCentMin(0.),
75fCentMax(100.),
76fNInputTracksMin(0),
77fNInputTracksMax(-1),
ea693273 78fAngStructCloseTracks(0),
da93bb11 79fCheckMethods(0),
75bf77e3 80fJetEtaMin(-.5),
81fJetEtaMax(.5),
82fJetPtMin(20.),
83fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
84fJetPtFractionMin(0.5),
85fNMatchJets(4),
86fMatchMaxDist(0.8),
87fKeepJets(kFALSE),
75bf77e3 88fkNbranches(2),
89fkEvtClasses(12),
90fOutputList(0x0),
91fbEvent(kTRUE),
92fHistEvtSelection(0x0),
da93bb11 93fhnDeltaR(0x0),
94fhnSumBkg(0x0),
95fhnJetCoreMethod3(0x0),
75bf77e3 96fh2JetCoreMethod1C10(0x0),
97fh2JetCoreMethod2C10(0x0),
ea693273 98fh2JetCoreMethod1C20(0x0),
99fh2JetCoreMethod2C20(0x0),
75bf77e3 100fh2JetCoreMethod1C30(0x0),
101fh2JetCoreMethod2C30(0x0),
75bf77e3 102fh2JetCoreMethod1C60(0x0),
103fh2JetCoreMethod2C60(0x0),
ea693273 104fh2AngStructpt1C10(0x0),
105fh2AngStructpt2C10(0x0),
106fh2AngStructpt3C10(0x0),
107fh2AngStructpt4C10(0x0),
108fh2AngStructpt1C20(0x0),
109fh2AngStructpt2C20(0x0),
110fh2AngStructpt3C20(0x0),
111fh2AngStructpt4C20(0x0),
112fh2AngStructpt1C30(0x0),
113fh2AngStructpt2C30(0x0),
114fh2AngStructpt3C30(0x0),
115fh2AngStructpt4C30(0x0),
116fh2AngStructpt1C60(0x0),
117fh2AngStructpt2C60(0x0),
118fh2AngStructpt3C60(0x0),
da93bb11 119fh2AngStructpt4C60(0x0),
120fh3spectriggered(0x0),
121fh3specbiased(0x0),
122fh3specleadsublead(0x0)
123
75bf77e3 124{
125 // default Constructor
126
127 fJetBranchName[0] = "";
128 fJetBranchName[1] = "";
129
130 fListJets[0] = new TList;
131 fListJets[1] = new TList;
132}
133
134AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
135AliAnalysisTaskSE(name),
136fESD(0x0),
137fAOD(0x0),
ea693273 138fAODExtension(0x0),
75bf77e3 139fBackgroundBranch(""),
ea693273 140fNonStdFile(""),
75bf77e3 141fIsPbPb(kTRUE),
142fOfflineTrgMask(AliVEvent::kAny),
143fMinContribVtx(1),
a9e585a7 144fVtxZMin(-10.),
145fVtxZMax(10.),
75bf77e3 146fEvtClassMin(0),
147fEvtClassMax(4),
8b47ec90 148fFilterMask(0),
ea693273 149fRadioFrac(0.2),
150fMinDist(0.1),
75bf77e3 151fCentMin(0.),
152fCentMax(100.),
153fNInputTracksMin(0),
154fNInputTracksMax(-1),
ea693273 155fAngStructCloseTracks(0),
da93bb11 156fCheckMethods(0),
75bf77e3 157fJetEtaMin(-.5),
158fJetEtaMax(.5),
159fJetPtMin(20.),
160fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
161fJetPtFractionMin(0.5),
162fNMatchJets(4),
163fMatchMaxDist(0.8),
164fKeepJets(kFALSE),
75bf77e3 165fkNbranches(2),
166fkEvtClasses(12),
167fOutputList(0x0),
168fbEvent(kTRUE),
169fHistEvtSelection(0x0),
da93bb11 170fhnDeltaR(0x0),
171fhnSumBkg(0x0),
172fhnJetCoreMethod3(0x0),
75bf77e3 173fh2JetCoreMethod1C10(0x0),
174fh2JetCoreMethod2C10(0x0),
ea693273 175fh2JetCoreMethod1C20(0x0),
176fh2JetCoreMethod2C20(0x0),
75bf77e3 177fh2JetCoreMethod1C30(0x0),
178fh2JetCoreMethod2C30(0x0),
75bf77e3 179fh2JetCoreMethod1C60(0x0),
180fh2JetCoreMethod2C60(0x0),
ea693273 181fh2AngStructpt1C10(0x0),
182fh2AngStructpt2C10(0x0),
183fh2AngStructpt3C10(0x0),
184fh2AngStructpt4C10(0x0),
185fh2AngStructpt1C20(0x0),
186fh2AngStructpt2C20(0x0),
187fh2AngStructpt3C20(0x0),
188fh2AngStructpt4C20(0x0),
189fh2AngStructpt1C30(0x0),
190fh2AngStructpt2C30(0x0),
191fh2AngStructpt3C30(0x0),
192fh2AngStructpt4C30(0x0),
193fh2AngStructpt1C60(0x0),
194fh2AngStructpt2C60(0x0),
195fh2AngStructpt3C60(0x0),
da93bb11 196fh2AngStructpt4C60(0x0),
197fh3spectriggered(0x0),
198fh3specbiased(0x0),
199fh3specleadsublead(0x0)
75bf77e3 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
212AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
213{
214 delete fListJets[0];
215 delete fListJets[1];
216}
217
218void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
219{
220 fJetBranchName[0] = branch1;
221 fJetBranchName[1] = branch2;
222}
223
224void 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
234void 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
da93bb11 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);
ea693273 259
da93bb11 260 if(fCheckMethods){
261 entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 ;
262 fhnJetCoreMethod3 = NewTHnSparseF("fhnEvent", entries);
75bf77e3 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);
ea693273 266 fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
267 fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
75bf77e3 268 fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
269 fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
75bf77e3 270 fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
da93bb11 271 fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
272
ea693273 273
da93bb11 274 if(fAngStructCloseTracks>0){
a9e585a7 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.); }
da93bb11 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);
75bf77e3 294
295 fOutputList->Add(fHistEvtSelection);
a9e585a7 296
da93bb11 297 fOutputList->Add(fhnDeltaR);
298
299 //fOutputList->Add(fhnSumBkg);
300
301
302
303 if(fCheckMethods){
304 fOutputList->Add(fhnJetCoreMethod3);
75bf77e3 305 fOutputList->Add(fh2JetCoreMethod1C10);
306 fOutputList->Add(fh2JetCoreMethod2C10);
ea693273 307 fOutputList->Add(fh2JetCoreMethod1C20);
308 fOutputList->Add(fh2JetCoreMethod2C20);
75bf77e3 309 fOutputList->Add(fh2JetCoreMethod1C30);
310 fOutputList->Add(fh2JetCoreMethod2C30);
75bf77e3 311 fOutputList->Add(fh2JetCoreMethod1C60);
da93bb11 312 fOutputList->Add(fh2JetCoreMethod2C60);}
a9e585a7 313
ea693273 314
75bf77e3 315
a9e585a7 316
317
318 if(fAngStructCloseTracks>0){
ea693273 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);
a9e585a7 334 fOutputList->Add(fh2AngStructpt4C60);}
da93bb11 335 fOutputList->Add(fh3spectriggered);
336 fOutputList->Add(fh3specbiased);
337 fOutputList->Add(fh3specleadsublead);
75bf77e3 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 }
ea693273 346 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
347 if (hn){
348 hn->Sumw2();
349 }
75bf77e3 350 }
351 TH1::AddDirectory(oldStatus);
352
353 PostData(1, fOutputList);
354}
355
356void 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 }
ea693273 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
75bf77e3 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());
da93bb11 391 cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
75bf77e3 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;
da93bb11 422 Double_t centValue = 0.;
75bf77e3 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
568f8fa2 433 fHistEvtSelection->Fill(0);
434 // accepted events
75bf77e3 435 // -- end event selection --
ea693273 436
75bf77e3 437 // get background
438 AliAODJetEventBackground* externalBackground = 0;
ea693273 439 if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
75bf77e3 440 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
441 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
442 }
ea693273 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
75bf77e3 448 Float_t rho = 0;
449 if(externalBackground)rho = externalBackground->GetBackground(0);
450
451
452 // fetch jets
453 TClonesArray *aodJets[2];
ea693273 454 aodJets[0]=0;
455 if(fAOD&&!aodJets[0]){
a9e585a7 456 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data()));
ea693273 457 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); }
458 if(fAODExtension && !aodJets[0]){
a9e585a7 459 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
ea693273 460 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
461
a9e585a7 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
ea693273 468 TList ParticleList;
469 Int_t nT = GetListOfTracks(&ParticleList);
470 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
75bf77e3 471 fListJets[iJetType]->Clear();
472 if (!aodJets[iJetType]) continue;
473
474 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
475
a9e585a7 476
75bf77e3 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);
a9e585a7 480 if(iJetType==0){
481 ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();
482 }}}
75bf77e3 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;
da93bb11 491 //Double_t distr=0.;
75bf77e3 492 Double_t phismall=0.;
a9e585a7 493 Int_t indexlead=-1;
494 Int_t indexsublead=-1;
495 Int_t indexstop=-1;
496
8b47ec90 497 if(fListJets[0]->GetEntries()>0) TMath::Sort(fListJets[0]->GetEntries(),ptsub,inord);
a9e585a7 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
da93bb11 514
ea693273 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.};
da93bb11 523
524
ea693273 525
75bf77e3 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();
ea693273 533 Double_t ptcorr=ptbig-rho*areabig;
534 if(ptcorr<=0) continue;
75bf77e3 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;
da93bb11 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
75bf77e3 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();
ea693273 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;
75bf77e3 583 index2=j;}
ea693273 584 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
585 index1=j;}} //en of loop over R=0.2 jets
75bf77e3 586 //method1:most concentric jet=core
587 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
ea693273 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); }
75bf77e3 592 //method2:hardest contained jet=core
593 if(index2!=-1){
594 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
ea693273 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);
da93bb11 598 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
ea693273 599
da93bb11 600
ea693273 601
da93bb11 602 for(int it = 0;it<nT;++it){
ea693273 603 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
da93bb11 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 }
568f8fa2 612 //end of track loop
da93bb11 613
614 //fhnSumBkg->Fill(centValue,ptcorr,bkg/jetbig->Pt(),partback->Pt(),flaglead);
615
616
ea693273 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
a9e585a7 647 if(fAngStructCloseTracks==2){
ea693273 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;
da93bb11 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));
ea693273 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;}}}}}
75bf77e3 680
75bf77e3 681
682
ea693273 683
684
685
686
687
688
689
690
691
692
693
694
75bf77e3 695 }
ea693273 696
697
698 //end loop over R=0.4 jets
da93bb11 699 if(fAngStructCloseTracks>0){
ea693273 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]);
a9e585a7 721 if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
ea693273 722
723
724
725
726
75bf77e3 727
728
729 PostData(1, fOutputList);
da93bb11 730}
75bf77e3 731
732void 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
ea693273 741
742
743
744
745
746
747
748
749
750
751Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
752
753 Int_t iCount = 0;
754
8b47ec90 755
ea693273 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
da93bb11 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
ea693273 803 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
804
805 Int_t iCount = 0;
806
8b47ec90 807
ea693273 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
75bf77e3 835Int_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
ea693273 872Double_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
75bf77e3 884
885
da93bb11 886THnSparse* 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
921void 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