]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetAntenna.cxx
CommitLineData
e59c88bd 1// ******************************************
2// This task computes several jet observables like
3// the fraction of energy in inner and outer coronnas,
4// jet-track correlations,triggered jet shapes and
5// correlation strength distribution of particles inside jets.
6// Author: lcunquei@cern.ch
7// *******************************************
8
9
10/**************************************************************************
11 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
12 * *
13 * Author: The ALICE Off-line Project. *
14 * Contributors are mentioned in the code where appropriate. *
15 * *
16 * Permission to use, copy, modify and distribute this software and its *
17 * documentation strictly for non-commercial purposes is hereby granted *
18 * without fee, provided that the above copyright notice appears in all *
19 * copies and that both the copyright notice and this permission notice *
20 * appear in the supporting documentation. The authors make no claims *
21 * about the suitability of this software for any purpose. It is *
22 * provided "as is" without express or implied warranty. *
23 **************************************************************************/
24
25
26#include "TChain.h"
27#include "TTree.h"
28#include "TMath.h"
29#include "TH1F.h"
30#include "TH2F.h"
31#include "TH3F.h"
32#include "THnSparse.h"
33#include "TMatrixD.h"
34#include "TMatrixDSym.h"
35#include "TMatrixDSymEigen.h"
36#include "TCanvas.h"
37#include "TRandom3.h"
38#include "AliLog.h"
39
40#include "AliAnalysisTask.h"
41#include "AliAnalysisManager.h"
42
43#include "AliVEvent.h"
44#include "AliESDEvent.h"
45#include "AliESDInputHandler.h"
46#include "AliCentrality.h"
47#include "AliAnalysisHelperJetTasks.h"
48#include "AliInputEventHandler.h"
49#include "AliAODJetEventBackground.h"
50#include "AliAODMCParticle.h"
51#include "AliAnalysisTaskFastEmbedding.h"
52#include "AliAODEvent.h"
53#include "AliAODHandler.h"
54#include "AliAODJet.h"
55
56#include "AliAnalysisTaskJetAntenna.h"
57
58using std::cout;
59using std::endl;
60
61ClassImp(AliAnalysisTaskJetAntenna)
62
63AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna() :
64AliAnalysisTaskSE(),
65fESD(0x0),
66fAODIn(0x0),
67fAODOut(0x0),
68fAODExtension(0x0),
69fBackgroundBranch(""),
70fNonStdFile(""),
71fIsPbPb(kTRUE),
72fOfflineTrgMask(AliVEvent::kAny),
73fMinContribVtx(1),
74fVtxZMin(-10.),
75fVtxZMax(10.),
76fEvtClassMin(0),
77fEvtClassMax(4),
78fFilterMask(0),
79fFilterMaskBestPt(0),
80fFilterType(0),
81fCentMin(0.),
82fCentMax(100.),
83fNInputTracksMin(0),
84fNInputTracksMax(-1),
85fRequireITSRefit(0),
86fApplySharedClusterCut(0),
87fTrackTypeRec(kTrackUndef),
88fRPAngle(0),
89fNRPBins(50),
90fSemigoodCorrect(0),
54aaacb8 91fDoMatching(kFALSE),
e59c88bd 92fHolePos(4.71),
93fHoleWidth(0.2),
94fCutTM(0.15),
95fJetEtaMin(-.5),
96fJetEtaMax(.5),
97fNevents(0),
98fTindex(0),
99fJetPtMin(20.),
100fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
101fJetPtFractionMin(0.5),
102fNMatchJets(4),
103fMatchMaxDist(0.8),
104fKeepJets(kFALSE),
105fkNbranches(2),
106fkEvtClasses(12),
107fOutputList(0x0),
108fHistEvtSelection(0x0),
109fh1JetEntries(0x0),
110fh2Circularity(0x0),
111fhnJetTM(0x0)
112{
113 // default Constructor
114
115 fJetBranchName[0] = "";
116 fJetBranchName[1] = "";
117
118 fListJets[0] = new TList;
119 fListJets[1] = new TList;
120}
121
122AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna(const char *name) :
123AliAnalysisTaskSE(name),
124fESD(0x0),
125fAODIn(0x0),
126fAODOut(0x0),
127fAODExtension(0x0),
128fBackgroundBranch(""),
129fNonStdFile(""),
130fIsPbPb(kTRUE),
131fOfflineTrgMask(AliVEvent::kAny),
132fMinContribVtx(1),
133fVtxZMin(-10.),
134fVtxZMax(10.),
135fEvtClassMin(0),
136fEvtClassMax(4),
137fFilterMask(0),
138fFilterMaskBestPt(0),
139fFilterType(0),
140fCentMin(0.),
141fCentMax(100.),
142fNInputTracksMin(0),
143fNInputTracksMax(-1),
144fRequireITSRefit(0),
145fApplySharedClusterCut(0),
146fTrackTypeRec(kTrackUndef),
147fRPAngle(0),
148fNRPBins(50),
149fSemigoodCorrect(0),
54aaacb8 150fDoMatching(kFALSE),
e59c88bd 151fHolePos(4.71),
152fHoleWidth(0.2),
153fCutTM(0.15),
154fJetEtaMin(-.5),
155fJetEtaMax(.5),
156fNevents(0),
157fTindex(0),
158fJetPtMin(20.),
159fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
160fJetPtFractionMin(0.5),
161fNMatchJets(4),
162fMatchMaxDist(0.8),
163fKeepJets(kFALSE),
164fkNbranches(2),
165fkEvtClasses(12),
166fOutputList(0x0),
167fHistEvtSelection(0x0),
168fh1JetEntries(0x0),
169fh2Circularity(0x0),
170fhnJetTM(0x0)
171 {
172 // Constructor
173
174
175 fJetBranchName[0] = "";
176 fJetBranchName[1] = "";
177
178 fListJets[0] = new TList;
179 fListJets[1] = new TList;
180
181 DefineOutput(1, TList::Class());
182}
183
184AliAnalysisTaskJetAntenna::~AliAnalysisTaskJetAntenna()
185{
186 delete fListJets[0];
187 delete fListJets[1];
188}
189
190void AliAnalysisTaskJetAntenna::SetBranchNames(const TString &branch1, const TString &branch2)
191{
192 fJetBranchName[0] = branch1;
193 fJetBranchName[1] = branch2;
194}
195
196void AliAnalysisTaskJetAntenna::Init()
197{
198
199 // check for jet branches
200 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
201 AliError("Jet branch name not set.");
202 }
203
204}
205
206void AliAnalysisTaskJetAntenna::UserCreateOutputObjects()
207{
208 // Create histograms
209 // Called once
210 OpenFile(1);
211 if(!fOutputList) fOutputList = new TList;
212 fOutputList->SetOwner(kTRUE);
213
214 Bool_t oldStatus = TH1::AddDirectoryStatus();
215 TH1::AddDirectory(kFALSE);
216
217 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
218 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
219 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
220 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
221 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
222 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
223 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
224 fOutputList->Add(fHistEvtSelection);
225 fh1JetEntries=new TH1F("JetEntries","",150,0,150);
226 fOutputList->Add(fh1JetEntries);
227 fh2Circularity=new TH2F("Circcularity","",10,0,1,150,0,150);
228 fOutputList->Add(fh2Circularity);
54aaacb8 229 Int_t nbinsJet[6]={15,30,9,36,10,20};
e59c88bd 230 Double_t binlowJet[6]= {0, 0, 0,-0.5*TMath::Pi(),0,0};
231 Double_t binupJet[6]= {1.5, 150,150,1.5*TMath::Pi(),1,200};
232 fhnJetTM = new THnSparseF("fhnJetTM", "fhnJetTM; dr;pt_jet;pt_track;phi;",6,nbinsJet,binlowJet,binupJet);
233 Double_t xPt3[10];
234 xPt3[0] = 0.;
235 for(Int_t i = 1;i<=9;i++){
236 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
237 else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
238 else xPt3[i] = xPt3[i-1] + 150.; // 18
239 }
240 fhnJetTM->SetBinEdges(2,xPt3);
241 fOutputList->Add(fhnJetTM);
242
243 // =========== Switch on Sumw2 for all histos ===========
244 for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
245 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
246 if (h1){
247 h1->Sumw2();
248 continue;
249 }
250 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
251 if (hn){
252 hn->Sumw2();
253 }
254 }
255 TH1::AddDirectory(oldStatus);
256
257 PostData(1, fOutputList);
258}
259
260void AliAnalysisTaskJetAntenna::UserExec(Option_t *)
261{
262
263
264 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
265 AliError("Jet branch name not set.");
266 return;
267 }
268
269 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
270 if (!fESD) {
271 AliError("ESD not available");
272 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
273 }
274 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
275
276 static AliAODEvent* aod = 0;
277 // take all other information from the aod we take the tracks from
278 if(!aod){
279 if(!fESD)aod = fAODIn;
280 else aod = fAODOut;}
281
282
283
284 if(fNonStdFile.Length()!=0){
285 // case that we have an AOD extension we need can fetch the jets from the extended output
286 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
287 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
288 if(!fAODExtension){
289 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
290 }
291 }
292
293
294 // -- event selection --
295 fHistEvtSelection->Fill(1); // number of events before event selection
296
297
298 Bool_t selected=kTRUE;
299 selected = AliAnalysisHelperJetTasks::Selected();
300 if(!selected){
301 // no selection by the service task, we continue
302 PostData(1,fOutputList);
303 return;}
304
305
306
307 // physics selection: this is now redundant, all should appear as accepted after service task selection
308 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
309 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
310 std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
311 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
312 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
313 fHistEvtSelection->Fill(2);
314 PostData(1, fOutputList);
315 return;
316 }
317
318
319
320 // vertex selection
321 if(!aod){
322 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
323 fHistEvtSelection->Fill(3);
324 PostData(1, fOutputList);
325 }
326 AliAODVertex* primVtx = aod->GetPrimaryVertex();
327
328 if(!primVtx){
329 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
330 fHistEvtSelection->Fill(3);
331 PostData(1, fOutputList);
332 return;
333 }
334
335 Int_t nTracksPrim = primVtx->GetNContributors();
336 if ((nTracksPrim < fMinContribVtx) ||
337 (primVtx->GetZ() < fVtxZMin) ||
338 (primVtx->GetZ() > fVtxZMax) ){
339 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
340 fHistEvtSelection->Fill(3);
341 PostData(1, fOutputList);
342 return;
343 }
344
345
346
347 // centrality selection
348 AliCentrality *cent = 0x0;
349 Double_t centValue = 0.;
350 if(fIsPbPb){
351 if(fESD) {cent = fESD->GetCentrality();
352 if(cent) centValue = cent->GetCentralityPercentile("V0M");}
353 else centValue=aod->GetHeader()->GetCentrality();
354
355 if(fDebug) printf("centrality: %f\n", centValue);
356 if (centValue < fCentMin || centValue > fCentMax){
357 fHistEvtSelection->Fill(4);
358 PostData(1, fOutputList);
359 return;
360 }}
361
362
363 fHistEvtSelection->Fill(0);
364 // accepted events
365 // -- end event selection --
366
367 // get background
368 AliAODJetEventBackground* externalBackground = 0;
369 if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
370 externalBackground = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
371 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
372 }
373 if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
374 externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
375 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
376 }
377
378 if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
379 externalBackground = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
380 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
381 }
382
383 Float_t rho = 0;
384
385
386 if(externalBackground)rho = externalBackground->GetBackground(0);
387
388 // fetch jets
389 TClonesArray *aodJets[2];
390 aodJets[0]=0;
391 if(fAODOut&&!aodJets[0]){
392 aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
393 aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));
394 }
395 if(fAODExtension && !aodJets[0]){
396 aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
397 aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
398 }
399 if(fAODIn&&!aodJets[0]){
400 aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
401 aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));
402 }
403
404
405
406 Int_t nT=0;
407 TList ParticleList;
54aaacb8 408 if(!fDoMatching) nT = GetListOfTracks(&ParticleList);
409 if(fDoMatching)nT=GetListOfTracksExtra(&ParticleList);
e59c88bd 410 if(nT<0){
411 PostData(1, fOutputList);
412 return;
413 }
414
415 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
416 fListJets[iJetType]->Clear();
417 if (!aodJets[iJetType]) continue;
418 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
419 for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
420 AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
421 if (jet) fListJets[iJetType]->Add(jet);
422 }
423 }
424
54aaacb8 425// jet matching
426 static TArrayI aMatchIndex(fListJets[0]->GetEntries());
427 static TArrayF aPtFraction(fListJets[0]->GetEntries());
428 if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
429 if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
e59c88bd 430
54aaacb8 431 if(fDoMatching){
432
e59c88bd 433
54aaacb8 434 // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
435 AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
436 fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
437 aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);}
438
439
440
441
442
443
444 // loop over matched jets
445 Int_t ir = -1; // index of matched reconstruced jet
446 Float_t fraction = -1.;
447 AliAODJet *jetmatched =0x0;
448
449
450
451
452
453 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
e59c88bd 454
455 Double_t etabig=0;
456 Double_t ptbig=0;
457 Double_t areabig=0;
458 Double_t phibig=0.;
459 Double_t pxbig,pybig,pzbig;
460
461
462 AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
463 etabig = jetbig->Eta();
464 phibig = jetbig->Phi();
465 ptbig = jetbig->Pt();
466 if(ptbig==0) continue;
467 areabig = jetbig->EffectiveAreaCharged();
54aaacb8 468 if(!fDoMatching) ptbig=ptbig-rho*areabig;
e59c88bd 469 if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
470
471 if(fSemigoodCorrect){
472 if((phibig>fHolePos-fHoleWidth) && (phibig<fHolePos+fHoleWidth)) continue;
473 }
e59c88bd 474 pxbig=jetbig->Px();
475 pybig=jetbig->Py();
476 pzbig=jetbig->Pz();
54aaacb8 477
478 if(fDoMatching){
479 Bool_t jetAccepted=kTRUE;
480 ir = aMatchIndex[i];
481 if(ir>=0) jetmatched = (AliAODJet*)(fListJets[1]->At(ir));
482 else continue;
483 fraction = aPtFraction[i];
484 // minimum fraction required
485
486 if(fraction<fJetPtFractionMin) jetAccepted = kFALSE;
487 if(jetAccepted){
488 // jet acceptance + minimum pT check
489 if(jetmatched->Eta()>fJetEtaMax || jetmatched->Eta()<fJetEtaMin) jetAccepted = kFALSE;
490
491 }
492 if(!jetAccepted) continue;
493 pxbig=jetmatched->Px();
494 pybig=jetmatched->Py();
495 pzbig=jetmatched->Pz();
496 ptbig=jetmatched->Pt()-rho*jetmatched->EffectiveAreaCharged();}
497
498
499
500
e59c88bd 501 TVector3 ppJ1(pxbig, pybig, pzbig);
502 TVector3 ppJ3(- pxbig * pzbig, - pybig * pzbig, pxbig * pxbig + pybig * pybig);
503 ppJ3.SetMag(1.);
504 TVector3 ppJ2(-pybig, pxbig, 0);
505 ppJ2.SetMag(1.);
506
54aaacb8 507
508
509
e59c88bd 510 Float_t mxx = 0.;
511 Float_t myy = 0.;
512 Float_t mxy = 0.;
513 Int_t nc = 0;
514 Float_t sump2 = 0.;
515
54aaacb8 516 for(int it = 0;it<ParticleList.GetEntries();++it){
e59c88bd 517 AliVParticle *track = (AliVParticle*)ParticleList.At(it);
518 TVector3 pp(track->Px(), track->Py(), track->Pz());
519 Float_t phi = track->Phi();
520 Float_t eta = track->Eta();
521 Float_t pt = track->Pt();
522 Float_t deta = eta - etabig;
523 Float_t dphi = RelativePhi(phi,phibig);
524 if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
525 Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
526 if (r < 0.4 && pt>fCutTM) {
527 //longitudinal and perpendicular component of the track pT in the
528 //local frame
529 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
530 TVector3 pPerp = pp - pLong;
531 //projection onto the two perpendicular vectors defined above
532 Float_t ppjX = pPerp.Dot(ppJ2);
533 Float_t ppjY = pPerp.Dot(ppJ3);
534 Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
535 //components of the 2D symmetrical sphericity matrix
536 mxx += (ppjX * ppjX / ppjT);
537 myy += (ppjY * ppjY / ppjT);
538 mxy += (ppjX * ppjY / ppjT);
539 nc++;
540 sump2 += ppjT;}
541 // max pt
542 if(nc<2) continue;
543
544 } // 1st Track Loop
545
546
547 // Sphericity Matrix
548 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
549 TMatrixDSym m0(2,ele);
550 // Find eigenvectors
551 TMatrixDSymEigen m(m0);
552 TVectorD eval(2);
553 TMatrixD evecm = m.GetEigenVectors();
554 eval = m.GetEigenValues();
555 // Largest eigenvector
556 Int_t jev = 0;
557 if (eval[0] < eval[1]) jev = 1;
558 TVectorD evec0(2);
559 // Principle axis
560 evec0 = TMatrixDColumn(evecm, jev);
561 TVector2 evec(evec0[0], evec0[1]);
562 Float_t circ=0;
563 if(jev==1) circ=2*eval[0];
564 if(jev==0) circ=2*eval[1];
565 fh2Circularity->Fill(circ,ptbig);
566 fh1JetEntries->Fill(ptbig);
54aaacb8 567 for (Int_t ip = 0; ip < ParticleList.GetEntries(); ip++) {
e59c88bd 568 AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
569 TVector3 pp(track->Px(), track->Py(), track->Pz());
570 Float_t phi = track->Phi();
571 Float_t eta = track->Eta();
572 Float_t pt = track->Pt();
573
574 Float_t deta = eta - etabig;
575 Float_t dphi = RelativePhi(phi,phibig);
576 if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
577
578 Float_t dRR = TMath::Sqrt(dphi * dphi + deta * deta);
579 TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
580 TVector3 pPerp = pp - pLong;
581 Float_t ppjX = pPerp.Dot(ppJ2);
582 Float_t ppjY = pPerp.Dot(ppJ3);
583 TVector2 vr(ppjX, ppjY) ;
584 //and this is the angle between the particle and the TM axis.
585 float phistr=evec.Phi()-vr.Phi();
586 if(phistr>2*TMath::Pi()) phistr -= 2*TMath::Pi();
587 if(phistr<-2*TMath::Pi()) phistr += 2*TMath::Pi();
588 if(phistr<-0.5*TMath::Pi()) phistr += 2*TMath::Pi();
589 if(phistr>1.5*TMath::Pi()) phistr -= 2*TMath::Pi();
590
4b3cfde0 591 double jetEntries[6] = {dRR,ptbig,pt,phistr,circ,static_cast<double>(nc)};
e59c88bd 592 fhnJetTM->Fill(jetEntries);
593
594 } // 2nd Track loop
595 }//jet loop
596
597 PostData(1, fOutputList);
598}
599
600void AliAnalysisTaskJetAntenna::Terminate(const Option_t *)
601{
602 // Draw result to the screen
603 // Called once at the end of the query
604
605 if (!GetOutputData(1))
606 return;
607}
608
609
610Int_t AliAnalysisTaskJetAntenna::GetListOfTracks(TList *list){
611
612 Int_t iCount = 0;
613 AliAODEvent *aod = 0;
614
615 if(!fESD)aod = fAODIn;
616 else aod = fAODOut;
e59c88bd 617 if(!aod)return 0;
618
e59c88bd 619
54aaacb8 620 for(int it = 0;it < aod->GetNumberOfTracks();++it){
e59c88bd 621 AliAODTrack *tr = aod->GetTrack(it);
622 Bool_t bGood = false;
623 if(fFilterType == 0)bGood = true;
624 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
625 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
626 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
627 if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
628 if(bGood==false) continue;
629 if (fApplySharedClusterCut) {
630 Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
631 if (frac > 0.4) continue;
632 }
633 if(TMath::Abs(tr->Eta())>0.9)continue;
634 if(tr->Pt()<0.15)continue;
635 list->Add(tr);
636 iCount++;
e59c88bd 637 }
54aaacb8 638 return iCount;
e59c88bd 639}
640
641
54aaacb8 642
643
644Int_t AliAnalysisTaskJetAntenna::GetListOfTracksExtra(TList *list){
645
646 Int_t iCount = 0;
647 AliAODEvent *aod = 0;
648
649 if(!fESD)aod = fAODIn;
650 else aod = fAODOut;
651 if(!aod)return 0;
652
653 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
654 if(!aodExtraTracks)return iCount;
655 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
656 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
657 if (!track) continue;
658 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
659 if(!trackAOD)continue;
660 Bool_t bGood = false;
661 if(fFilterType == 0)bGood = true;
662 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
663 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
664 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
665 if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
666 if (fApplySharedClusterCut) {
667 Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
668 if (frac > 0.4) continue;
669 }
670
671
672 if(TMath::Abs(trackAOD->Eta())>0.9) continue;
673 if(trackAOD->Pt()<0.15) continue;
674 if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
675 list->Add(trackAOD);
676 iCount++;
677 }
678
679 return iCount;
680}
681
e59c88bd 682Double_t AliAnalysisTaskJetAntenna::RelativePhi(Double_t mphi,Double_t vphi){
683
684 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
685 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
686 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
687 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
688 double dphi = mphi-vphi;
689 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
690 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
691 return dphi;//dphi in [-Pi, Pi]
692}
693
694Int_t AliAnalysisTaskJetAntenna::GetPhiBin(Double_t phi)
695{
696 Int_t phibin=-1;
697 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
698 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
699 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
700 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
701 return phibin;
702}