git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity
[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
bc684956 400 if(!fAOD){
401 PostData(1, fOutputList);
402 if(fDebug) Printf(" No AOD found ... ");
403 return;
404 }
75bf77e3 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;
da93bb11 427 Double_t centValue = 0.;
75bf77e3 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
568f8fa2 438 fHistEvtSelection->Fill(0);
439 // accepted events
75bf77e3 440 // -- end event selection --
ea693273 441
75bf77e3 442 // get background
443 AliAODJetEventBackground* externalBackground = 0;
ea693273 444 if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
75bf77e3 445 externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
446 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
447 }
ea693273 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
75bf77e3 453 Float_t rho = 0;
454 if(externalBackground)rho = externalBackground->GetBackground(0);
455
456
457 // fetch jets
458 TClonesArray *aodJets[2];
ea693273 459 aodJets[0]=0;
460 if(fAOD&&!aodJets[0]){
a9e585a7 461 aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data()));
ea693273 462 aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); }
463 if(fAODExtension && !aodJets[0]){
a9e585a7 464 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
ea693273 465 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data())); }
466
a9e585a7 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
ea693273 473 TList ParticleList;
474 Int_t nT = GetListOfTracks(&ParticleList);
475 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
75bf77e3 476 fListJets[iJetType]->Clear();
477 if (!aodJets[iJetType]) continue;
75bf77e3 478 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
479
a9e585a7 480
75bf77e3 481 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
482 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
bc684956 483 if(!jet)continue;
484 fListJets[iJetType]->Add(jet);
a9e585a7 485 if(iJetType==0){
486 ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();
487 }}}
75bf77e3 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;
da93bb11 496 //Double_t distr=0.;
75bf77e3 497 Double_t phismall=0.;
a9e585a7 498 Int_t indexlead=-1;
499 Int_t indexsublead=-1;
500 Int_t indexstop=-1;
501
8b47ec90 502 if(fListJets[0]->GetEntries()>0) TMath::Sort(fListJets[0]->GetEntries(),ptsub,inord);
a9e585a7 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
da93bb11 519
ea693273 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.};
da93bb11 528
529
ea693273 530
75bf77e3 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();
ea693273 538 Double_t ptcorr=ptbig-rho*areabig;
539 if(ptcorr<=0) continue;
75bf77e3 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;
da93bb11 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
75bf77e3 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();
ea693273 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;
75bf77e3 588 index2=j;}
ea693273 589 if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
590 index1=j;}} //en of loop over R=0.2 jets
75bf77e3 591 //method1:most concentric jet=core
592 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
ea693273 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); }
75bf77e3 597 //method2:hardest contained jet=core
598 if(index2!=-1){
599 AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
ea693273 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);
da93bb11 603 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
ea693273 604
da93bb11 605
ea693273 606
da93bb11 607 for(int it = 0;it<nT;++it){
ea693273 608 AliVParticle *part = (AliVParticle*)ParticleList.At(it);
da93bb11 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 }
568f8fa2 617 //end of track loop
da93bb11 618
619 //fhnSumBkg->Fill(centValue,ptcorr,bkg/jetbig->Pt(),partback->Pt(),flaglead);
620
621
ea693273 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
a9e585a7 652 if(fAngStructCloseTracks==2){
ea693273 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;
da93bb11 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));
ea693273 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;}}}}}
75bf77e3 685
75bf77e3 686
687
ea693273 688
689
690
691
692
693
694
695
696
697
698
699
75bf77e3 700 }
ea693273 701
702
703 //end loop over R=0.4 jets
da93bb11 704 if(fAngStructCloseTracks>0){
ea693273 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]);
a9e585a7 726 if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
ea693273 727
728
729
730
731
75bf77e3 732
733
734 PostData(1, fOutputList);
da93bb11 735}
75bf77e3 736
737void 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
ea693273 746
747
748
749
750
751
752
753
754
755
756Int_t AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
757
758 Int_t iCount = 0;
759
8b47ec90 760
ea693273 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
da93bb11 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
ea693273 808 Int_t AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
809
810 Int_t iCount = 0;
811
8b47ec90 812
ea693273 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
75bf77e3 840Int_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
ea693273 877Double_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
75bf77e3 889
890
da93bb11 891THnSparse* 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
926void 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;
bc684956 988 break;
da93bb11 989
da93bb11 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