]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliAnalysisTaskJetCorePP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
additional fixes from Salvatore
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
CommitLineData
ad869500 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 "TH1D.h"
32#include "TH2F.h"
33#include "TH3F.h"
34#include "THnSparse.h"
35#include "TCanvas.h"
36
37#include "AliLog.h"
38
39#include "AliAnalysisTask.h"
40#include "AliAnalysisManager.h"
41
42#include "AliVEvent.h"
43#include "AliESDEvent.h"
44#include "AliESDInputHandler.h"
45#include "AliCentrality.h"
46#include "AliAnalysisHelperJetTasks.h"
47#include "AliInputEventHandler.h"
48#include "AliAODJetEventBackground.h"
49#include "AliAODMCParticle.h"
50//#include "AliAnalysisTaskFastEmbedding.h"
51#include "AliAODEvent.h"
52#include "AliAODHandler.h"
53#include "AliAODJet.h"
54
55#include "AliAnalysisTaskJetCorePP.h"
56
57using std::cout;
58using std::endl;
59
60ClassImp(AliAnalysisTaskJetCorePP)
61
62//Filip Krizek 1st March 2013
63
64//---------------------------------------------------------------------
65AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
66AliAnalysisTaskSE(),
67fESD(0x0),
68fAODIn(0x0),
69fAODOut(0x0),
70fAODExtension(0x0),
71fJetBranchName(""),
7fc3d134 72fListJets(0x0),
ad869500 73fNonStdFile(""),
74fSystem(0), //pp=0 pPb=1
75fJetParamR(0.4),
76fOfflineTrgMask(AliVEvent::kAny),
77fMinContribVtx(1),
78fVtxZMin(-10.0),
79fVtxZMax(10.0),
80fFilterMask(0),
81fCentMin(0.0),
82fCentMax(100.0),
83fJetEtaMin(-0.5),
84fJetEtaMax(0.5),
85fTriggerEtaCut(0.9),
86fTrackEtaCut(0.9),
87fTrackLowPtCut(0.15),
88fOutputList(0x0),
89fHistEvtSelection(0x0),
90fh2Ntriggers(0x0),
91fHJetSpec(0x0),
92fHJetDensity(0x0),
93fHJetDensityA4(0x0),
94fhJetPhi(0x0),
95fhTriggerPhi(0x0),
96fhJetEta(0x0),
97fhTriggerEta(0x0),
98fhVertexZ(0x0),
99fhVertexZAccept(0x0),
100fhContribVtx(0x0),
101fhContribVtxAccept(0x0),
102fhDphiTriggerJet(0x0),
103fhDphiTriggerJetAccept(0x0),
104fhCentrality(0x0),
105fhCentralityAccept(0x0),
7fc3d134 106fHJetPtRaw(0x0),
107fHLeadingJetPtRaw(0x0),
108fHDphiVsJetPtAll(0x0),
109fHRhoFastJetVsRhoCone(0x0),
110fkAcceptance(2.0*TMath::Pi()*1.8),
111fConeArea(TMath::Pi()*0.4*0.4)
ad869500 112{
113 // default Constructor
ad869500 114}
115
116//---------------------------------------------------------------------
117
118AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
119AliAnalysisTaskSE(name),
120fESD(0x0),
121fAODIn(0x0),
122fAODOut(0x0),
123fAODExtension(0x0),
124fJetBranchName(""),
7fc3d134 125fListJets(0x0),
ad869500 126fNonStdFile(""),
127fSystem(0), //pp=0 pPb=1
128fJetParamR(0.4),
129fOfflineTrgMask(AliVEvent::kAny),
130fMinContribVtx(1),
131fVtxZMin(-10.0),
132fVtxZMax(10.0),
133fFilterMask(0),
134fCentMin(0.0),
135fCentMax(100.0),
136fJetEtaMin(-0.5),
137fJetEtaMax(0.5),
138fTriggerEtaCut(0.9),
139fTrackEtaCut(0.9),
140fTrackLowPtCut(0.15),
141fOutputList(0x0),
142fHistEvtSelection(0x0),
143fh2Ntriggers(0x0),
144fHJetSpec(0x0),
145fHJetDensity(0x0),
146fHJetDensityA4(0x0),
147fhJetPhi(0x0),
148fhTriggerPhi(0x0),
149fhJetEta(0x0),
150fhTriggerEta(0x0),
151fhVertexZ(0x0),
152fhVertexZAccept(0x0),
153fhContribVtx(0x0),
154fhContribVtxAccept(0x0),
155fhDphiTriggerJet(0x0),
156fhDphiTriggerJetAccept(0x0),
157fhCentrality(0x0),
158fhCentralityAccept(0x0),
7fc3d134 159fHJetPtRaw(0x0),
160fHLeadingJetPtRaw(0x0),
161fHDphiVsJetPtAll(0x0),
162fHRhoFastJetVsRhoCone(0x0),
163fkAcceptance(2.0*TMath::Pi()*1.8),
164fConeArea(TMath::Pi()*0.4*0.4)
ad869500 165{
7fc3d134 166// Constructor
ad869500 167
168 DefineOutput(1, TList::Class());
169}
170
171//--------------------------------------------------------------
172AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
173AliAnalysisTaskSE(a.GetName()),
174fESD(a.fESD),
175fAODIn(a.fAODIn),
176fAODOut(a.fAODOut),
177fAODExtension(a.fAODExtension),
178fJetBranchName(a.fJetBranchName),
179fListJets(a.fListJets),
180fNonStdFile(a.fNonStdFile),
181fSystem(a.fSystem),
182fJetParamR(a.fJetParamR),
183fOfflineTrgMask(a.fOfflineTrgMask),
184fMinContribVtx(a.fMinContribVtx),
185fVtxZMin(a.fVtxZMin),
186fVtxZMax(a.fVtxZMax),
187fFilterMask(a.fFilterMask),
188fCentMin(a.fCentMin),
189fCentMax(a.fCentMax),
190fJetEtaMin(a.fJetEtaMin),
191fJetEtaMax(a.fJetEtaMax),
192fTriggerEtaCut(a.fTriggerEtaCut),
193fTrackEtaCut(a.fTrackEtaCut),
194fTrackLowPtCut(a.fTrackLowPtCut),
195fOutputList(a.fOutputList),
196fHistEvtSelection(a.fHistEvtSelection),
197fh2Ntriggers(a.fh2Ntriggers),
198fHJetSpec(a.fHJetSpec),
199fHJetDensity(a.fHJetDensity),
200fHJetDensityA4(a.fHJetDensityA4),
201fhJetPhi(a.fhJetPhi),
202fhTriggerPhi(a.fhTriggerPhi),
203fhJetEta(a.fhJetEta),
204fhTriggerEta(a.fhTriggerEta),
205fhVertexZ(a.fhVertexZ),
206fhVertexZAccept(a.fhVertexZAccept),
207fhContribVtx(a.fhContribVtx),
208fhContribVtxAccept(a.fhContribVtxAccept),
209fhDphiTriggerJet(a.fhDphiTriggerJet),
210fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
211fhCentrality(a.fhCentrality),
212fhCentralityAccept(a.fhCentralityAccept),
7fc3d134 213fHJetPtRaw(a.fHJetPtRaw),
214fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
215fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
216fHRhoFastJetVsRhoCone(a.fHRhoFastJetVsRhoCone),
217fkAcceptance(a.fkAcceptance),
218fConeArea(a.fConeArea)
ad869500 219{
220 //Copy Constructor
221}
222//--------------------------------------------------------------
223
224AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
225 // assignment operator
226 this->~AliAnalysisTaskJetCorePP();
227 new(this) AliAnalysisTaskJetCorePP(a);
228 return *this;
229}
230//--------------------------------------------------------------
231
232AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
233{
234 //Destructor
235 delete fListJets;
236 delete fOutputList; // ????
237}
238
239//--------------------------------------------------------------
240
241void AliAnalysisTaskJetCorePP::Init()
242{
243 // check for jet branches
244 if(!strlen(fJetBranchName.Data())){
245 AliError("Jet branch name not set.");
246 }
247
248}
249
250//--------------------------------------------------------------
251
252void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
253{
254
255
256 // Create histograms
257 // Called once
7fc3d134 258 fListJets = new TList();
259
ad869500 260 OpenFile(1);
261 if(!fOutputList) fOutputList = new TList();
262 fOutputList->SetOwner(kTRUE);
263
264 Bool_t oldStatus = TH1::AddDirectoryStatus();
265 TH1::AddDirectory(kFALSE);
266
267 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
268 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
269 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
270 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
271 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
272 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
273 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
274
275 fOutputList->Add(fHistEvtSelection);
276
277 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
278
279 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
280 nBinsCentrality,0.0,100.0,50,0.0,50.0);
7fc3d134 281 fOutputList->Add(fh2Ntriggers);
ad869500 282
283 //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
284 const Int_t dimSpec = 6;
7fc3d134 285 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 140, 50, 100, 100};
ad869500 286 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0,-80.0, 0.0, 0.0, 0.0};
7fc3d134 287 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0,200.0,50.0, 200, 50.0};
288 fHJetSpec = new THnSparseF("fHJetSpec",
289 "Recoil jet spectrum [cent,A,pTjet-rho*A,pTtrig,pTjet, rho*A]",
290 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
ad869500 291 fOutputList->Add(fHJetSpec);
292
293 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
294 //Jet number density histos [Trk Mult, jet density, pT trigger]
295 const Int_t dimJetDens = 3;
296 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
297 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
298 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0 };
299
300 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
301 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
302
303 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
304 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
305
ad869500 306 fOutputList->Add(fHJetDensity);
307 fOutputList->Add(fHJetDensityA4);
308
309
310 //inclusive azimuthal and pseudorapidity histograms
7fc3d134 311 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
312 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
313 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
314 25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
315 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
316 50,0, 100, 40,-0.9,0.9);
317 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
318 25, 0, 50, 40,-0.9,0.9);
319
ad869500 320 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
321 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
322 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
323 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
324 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
325 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
326 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
327 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
328
329 fOutputList->Add(fhJetPhi);
330 fOutputList->Add(fhTriggerPhi);
331 fOutputList->Add(fhJetEta);
332 fOutputList->Add(fhTriggerEta);
333 fOutputList->Add(fhVertexZ);
334 fOutputList->Add(fhVertexZAccept);
335 fOutputList->Add(fhContribVtx);
336 fOutputList->Add(fhContribVtxAccept);
337 fOutputList->Add(fhDphiTriggerJet);
338 fOutputList->Add(fhDphiTriggerJetAccept);
339 fOutputList->Add(fhCentrality);
340 fOutputList->Add(fhCentralityAccept);
341
7fc3d134 342 // raw spectra of INCLUSIVE jets
343 //Centrality, pTjet, pTjet - rho*A, A
344 const Int_t dimRaw = 4;
345 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 75, 100};
346 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, -50.0, 0.0};
347 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 100.0, 1.0};
348 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
349 "Incl. jet spectrum [cent,pTjet,pTjet-rho*A,A]",
350 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
351 fOutputList->Add(fHJetPtRaw);
352
353 // raw spectra of LEADING jets
354 //Centrality, pTjet, pTjet - rho*A, A
355 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
356 "Leading jet spectrum [cent,pTjet,pTjet-rho*A,A]",
357 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
358 fOutputList->Add(fHLeadingJetPtRaw);
359
360 // Dphi versus pT jet
361 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
362 const Int_t dimDp = 4;
363 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
364 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
365 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
366 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
367 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
368 dimDp,nBinsDp,lowBinDp,hiBinDp);
369 fOutputList->Add(fHDphiVsJetPtAll);
370
371 // Rho Fast Jet vs Rho Cone
372 //Centrality, Rho Fast Jet, Rho Perp Cone, pTjet
373 const Int_t dimRo = 4;
374 const Int_t nBinsRo[dimRo] = {nBinsCentrality, 100, 100, 50};
375 const Double_t lowBinRo[dimRo] = {0.0, 0.0, 0.0, 0.0};
376 const Double_t hiBinRo[dimRo] = {100.0, 100.0, 100.0, 100.0};
377 fHRhoFastJetVsRhoCone = new THnSparseF("fHRhoFastJetVsRhoCone",
378 "Rho FastJet vs rho PerpCone [cent,rhoFastJet,rhoCone,pTjet]",
379 dimRo,nBinsRo,lowBinRo,hiBinRo);
380 fOutputList->Add(fHRhoFastJetVsRhoCone);
381
382
383
ad869500 384
385 // =========== Switch on Sumw2 for all histos ===========
386 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
387 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
388 if(h1){
389 h1->Sumw2();
390 continue;
391 }
392 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
393 if(hn){
394 hn->Sumw2();
395 }
396 }
397 TH1::AddDirectory(oldStatus);
398
399 PostData(1, fOutputList);
400}
401
402//--------------------------------------------------------------------
403
404void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
405{
406
407 //Event loop
408
409 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
410 AliError("Cone radius is set to zero.");
411 return;
412 }
413 if(!strlen(fJetBranchName.Data())){
414 AliError("Jet branch name not set.");
415 return;
416 }
417
418 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
419 if(!fESD){
420 AliError("ESD not available");
421 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
422 }
423
424 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
425 AliAODEvent* aod = NULL;
426 // take all other information from the aod we take the tracks from
427 if(!aod){
428 if(!fESD) aod = fAODIn;
429 else aod = fAODOut;
430 }
431
432 if(fNonStdFile.Length()!=0){
433 // case that we have an AOD extension we can fetch the jets from the extended output
434 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
435 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
436 if(!fAODExtension){
437 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
438 }
439 }
440
441 // ----------------- event selection --------------------------
442 fHistEvtSelection->Fill(1); // number of events before event selection
443
444 // physics selection
445 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
446 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
447
448 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
449 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
450 fHistEvtSelection->Fill(2);
451 PostData(1, fOutputList);
452 return;
453 }
454
455 //check AOD pointer
456 if(!aod){
457 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
458 fHistEvtSelection->Fill(3);
459 PostData(1, fOutputList);
460 return;
461 }
462
463 // vertex selection
464 AliAODVertex* primVtx = aod->GetPrimaryVertex();
465
466 if(!primVtx){
467 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
468 fHistEvtSelection->Fill(3);
469 PostData(1, fOutputList);
470 return;
471 }
472
473 Int_t nTracksPrim = primVtx->GetNContributors();
474 Float_t vtxz = primVtx->GetZ();
475 //Input events
476 fhContribVtx->Fill(nTracksPrim);
477 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
478
479 if((nTracksPrim < fMinContribVtx) ||
480 (vtxz < fVtxZMin) ||
481 (vtxz > fVtxZMax)){
482 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
483 (char*)__FILE__,__LINE__,vtxz);
484 fHistEvtSelection->Fill(3);
485 PostData(1, fOutputList);
486 return;
487 }else{
488 //Accepted events
489 fhContribVtxAccept->Fill(nTracksPrim);
490 fhVertexZAccept->Fill(vtxz);
491 }
492 //FK// No event class selection imposed
493 // event class selection (from jet helper task)
494 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
495 //if(fDebug) Printf("Event class %d", eventClass);
496 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
497 // fHistEvtSelection->Fill(4);
498 // PostData(1, fOutputList);
499 // return;
500 //}
501
502 // centrality selection
503 AliCentrality *cent = 0x0;
504 Double_t centValue = 0.0;
505 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
506 if(fESD){
507 cent = fESD->GetCentrality();
508 if(cent) centValue = cent->GetCentralityPercentile("V0M");
509 }else{
510 centValue = aod->GetHeader()->GetCentrality();
511 }
512 if(fDebug) printf("centrality: %f\n", centValue);
513 //Input events
514 fhCentrality->Fill(centValue);
515
516 if(centValue < fCentMin || centValue > fCentMax){
517 fHistEvtSelection->Fill(4);
518 PostData(1, fOutputList);
519 return;
520 }else{
521 //Accepted events
522 fhCentralityAccept->Fill( centValue );
523 }
524 }
525
526 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
527
528 fHistEvtSelection->Fill(0); // accepted events
7fc3d134 529 fConeArea = TMath::Pi()*fJetParamR*fJetParamR;
ad869500 530 // ------------------- end event selection --------------------
531
532 // fetch jets
533 TClonesArray *aodJets = 0x0;
534
535 if(fAODOut && !aodJets){
536 aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
537 }
538 if(fAODExtension && !aodJets){
539 aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
540 }
541 if(fAODIn && !aodJets){
542 aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
543 }
544
545 // ------------- Hadron trigger --------------
546 TList particleList; //list of tracks
547 Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
548
7fc3d134 549 if(indexTrigg<0) return; // no trigger track found above 150 MeV/c
ad869500 550
551 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
552 if(!triggerHadron){
553 PostData(1, fOutputList);
554 return;
555 }
556
557
558 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
559
ad869500 560 //Trigger Diagnostics---------------------------------
7fc3d134 561 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
562 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
ad869500 563
564 //--------- Fill list of jets --------------
565 fListJets->Clear();
566 if(aodJets){
567 if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(),
568 aodJets->GetEntriesFast());
569 for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
570 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
571 if (jet) fListJets->Add(jet);
572 }
573 }
574
575 Double_t etaJet = 0.0;
576 Double_t pTJet = 0.0;
577 Double_t areaJet = 0.0;
578 Double_t phiJet = 0.0;
7fc3d134 579 Int_t injet4 = 0;
580 Int_t injet = 0;
581 Int_t indexLeadingJet = -1;
582 Double_t pTLeadingJet = -10.0;
583 Double_t pTcorrLeadingJet = -10.0;
584 Double_t areaLeadingJet = -10.0;
585 Double_t rhoFastJet = 0.0;
ad869500 586 //---------- jet loop ---------
587 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
588 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
589 if(!jet) continue;
590 etaJet = jet->Eta();
591 phiJet = jet->Phi();
592 pTJet = jet->Pt();
593 if(pTJet==0) continue;
594
595 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
596 areaJet = jet->EffectiveAreaCharged();
7fc3d134 597
598 Double_t fastJetbgpT = jet->ChargedBgEnergy();
599 if(areaJet>0) rhoFastJet = fastJetbgpT/areaJet;
600 else rhoFastJet = 0.0;
601
ad869500 602 //Jet Diagnostics---------------------------------
7fc3d134 603 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
604 fhJetEta->Fill(pTJet, etaJet);
ad869500 605 if(areaJet >= 0.07) injet++;
606 if(areaJet >= 0.4) injet4++;
607 //--------------------------------------------------
608
609 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
610
611 fhDphiTriggerJet->Fill(dphi); //Input
ad869500 612 //Background w.r.t. jet axis
613 Double_t pTBckWrtJet =
614 GetBackgroundInPerpCone(fJetParamR, phiJet, etaJet, &particleList);
615
7fc3d134 616 Double_t ratioOfAreas = areaJet/fConeArea;
ad869500 617 Double_t rhoA = ratioOfAreas*pTBckWrtJet; //bg activity in a cone of similar size
618 Double_t ptcorr = pTJet - rhoA; //Pt Jet UE subtr
619
7fc3d134 620 //search for leading jet
621 if(pTJet > pTLeadingJet){
622 indexLeadingJet = ij;
623 pTLeadingJet = pTJet;
624 pTcorrLeadingJet = ptcorr;
625 areaLeadingJet = areaJet;
626 }
627
628 // raw spectra of INCLUSIVE jets
629 //Centrality, pTjet, pTjet - rho*A, A
630 Double_t fillraw[] = { centValue,
631 pTJet,
632 ptcorr,
633 areaJet
634 };
635 fHJetPtRaw->Fill(fillraw);
636
637 //Dphi versus jet pT
638 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
639 Double_t filldp[] = { centValue,
640 dphi,
641 pTJet,
642 triggerHadron->Pt()
643 };
644 fHDphiVsJetPtAll->Fill(filldp);
645
646 // Rho Fast Jet vs Rho Cone
647 //Centrality, Rho Fast Jet, Rho Perp Cone, pTjet
648 Double_t fillro[] = { centValue,
649 rhoFastJet,
650 pTBckWrtJet/fConeArea,
651 pTJet
652 };
653 fHRhoFastJetVsRhoCone->Fill(fillro);
654
655 // Select back to back trigger - jet pairs
656 if(TMath::Abs((Double_t) dphi) < TMath::Pi()-0.6) continue;
657 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
658
659
ad869500 660 //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
661 Double_t fillspec[] = { centValue,
662 areaJet,
663 ptcorr,
664 triggerHadron->Pt(),
665 pTJet,
666 rhoA
667 };
668 fHJetSpec->Fill(fillspec);
669
670 }//jet loop
7fc3d134 671
672 if(indexLeadingJet > -1){
673 // raw spectra of LEADING jets
674 //Centrality, pTjet, pTjet - rho*A, A
675 Double_t fillleading[] = { centValue,
676 pTLeadingJet,
677 pTcorrLeadingJet,
678 areaLeadingJet
679 };
680 fHLeadingJetPtRaw->Fill(fillleading);
681 }
682
683
ad869500 684 //Fill Jet Density In the Event A>0.07
685 if(injet>0){
686 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
687 injet/fkAcceptance,
688 triggerHadron->Pt()
689 };
690 fHJetDensity->Fill(filldens);
691 }
692
693 //Fill Jet Density In the Event A>0.4
694 if(injet4>0){
695 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
696 injet4/fkAcceptance,
697 triggerHadron->Pt()
698 };
699 fHJetDensityA4->Fill(filldens4);
700 }
701
702 PostData(1, fOutputList);
703}
704
705//----------------------------------------------------------------------------
706void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
707{
708 // Draw result to the screen
709 // Called once at the end of the query
710
711 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
712 if(!GetOutputData(1)) return;
713}
714
715
716//----------------------------------------------------------------------------
717Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
718 //Fill the list of accepted tracks (passed track cut)
719 //return consecutive index of the hardest ch hadron in the list
720 Int_t iCount = 0;
721 AliAODEvent *aodevt = NULL;
722
723 if(!fESD) aodevt = fAODIn;
724 else aodevt = fAODOut;
725
726 if(!aodevt) return -1;
727
728 Int_t index = -1; //index of the highest particle in the list
729 Double_t ptmax = -10;
730
731 for(int it = 0; it < aodevt->GetNumberOfTracks(); it++){
732 AliAODTrack *tr = aodevt->GetTrack(it);
733
734 //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
735 if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
736 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
737 if(tr->Pt() < fTrackLowPtCut) continue;
738 list->Add(tr);
739 if(tr->Pt()>ptmax){
740 ptmax = tr->Pt();
741 index = iCount;
742 }
743 iCount++;
744 }
745
746 if(index>-1){ //check pseudorapidity cut on trigger
747 AliAODTrack *trigger = (AliAODTrack*) list->At(index);
748 if(trigger && TMath::Abs((Float_t) trigger->Eta())< fTriggerEtaCut){ return index;}
749 return -1;
750 }else{
7fc3d134 751 return -1;
ad869500 752 }
753}
754
755//----------------------------------------------------------------------------
756
757Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
758 //calculate sum of track pT in the cone perpendicular in phi to the jet
759 //jetR = cone radius
760 // jetPhi, jetEta = direction of the jet
761 Int_t numberOfTrks = trkList->GetEntries();
762 Double_t pTsum = 0.0;
763 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
764 for(Int_t it=0; it<numberOfTrks; it++){
765 AliVParticle *trk = (AliVParticle*) trkList->At(it);
766 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
767 Double_t deta = trk->Eta()-jetEta;
768
769 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
770 pTsum += trk->Pt();
771 }
772 }
773
774 return pTsum;
775}
776
777//----------------------------------------------------------------------------
778
779Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
780 //Get relative azimuthal angle of two particles -pi to pi
781 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
782 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
783
784 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
785 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
786
787 Double_t dphi = mphi - vphi;
788 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
789 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
790
791 return dphi;//dphi in [-Pi, Pi]
792}
793
794