Get Extra Tracks from Embedding
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskJetAntenna.cxx
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
58 using std::cout;
59 using std::endl;
60
61 ClassImp(AliAnalysisTaskJetAntenna)
62
63 AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna() :
64 AliAnalysisTaskSE(),
65 fESD(0x0),
66 fAODIn(0x0),
67 fAODOut(0x0),
68 fAODExtension(0x0),
69 fBackgroundBranch(""),
70 fNonStdFile(""),
71 fIsPbPb(kTRUE),
72 fOfflineTrgMask(AliVEvent::kAny),
73 fMinContribVtx(1),
74 fVtxZMin(-10.),
75 fVtxZMax(10.),
76 fEvtClassMin(0),
77 fEvtClassMax(4),
78 fFilterMask(0),
79 fFilterMaskBestPt(0),
80 fFilterType(0),
81 fCentMin(0.),
82 fCentMax(100.),
83 fNInputTracksMin(0),
84 fNInputTracksMax(-1),
85 fRequireITSRefit(0),
86 fApplySharedClusterCut(0),
87 fTrackTypeRec(kTrackUndef),
88 fRPAngle(0),
89 fNRPBins(50),
90 fSemigoodCorrect(0),
91 fDoMatching(kFALSE),
92 fHolePos(4.71),
93 fHoleWidth(0.2),
94 fCutTM(0.15),
95 fJetEtaMin(-.5),
96 fJetEtaMax(.5),
97 fNevents(0),
98 fTindex(0),
99 fJetPtMin(20.),
100 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
101 fJetPtFractionMin(0.5),
102 fNMatchJets(4),
103 fMatchMaxDist(0.8),
104 fKeepJets(kFALSE),
105 fkNbranches(2),
106 fkEvtClasses(12),
107 fOutputList(0x0),
108 fHistEvtSelection(0x0),
109 fh1JetEntries(0x0),
110 fh2Circularity(0x0),
111 fhnJetTM(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
122 AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna(const char *name) :
123 AliAnalysisTaskSE(name),
124 fESD(0x0),
125 fAODIn(0x0),
126 fAODOut(0x0),
127 fAODExtension(0x0),
128 fBackgroundBranch(""),
129 fNonStdFile(""),
130 fIsPbPb(kTRUE),
131 fOfflineTrgMask(AliVEvent::kAny),
132 fMinContribVtx(1),
133 fVtxZMin(-10.),
134 fVtxZMax(10.),
135 fEvtClassMin(0),
136 fEvtClassMax(4),
137 fFilterMask(0),
138 fFilterMaskBestPt(0),
139 fFilterType(0),
140 fCentMin(0.),
141 fCentMax(100.),
142 fNInputTracksMin(0),
143 fNInputTracksMax(-1),
144 fRequireITSRefit(0),
145 fApplySharedClusterCut(0),
146 fTrackTypeRec(kTrackUndef),
147 fRPAngle(0),
148 fNRPBins(50),
149 fSemigoodCorrect(0),
150 fDoMatching(kFALSE),
151 fHolePos(4.71),
152 fHoleWidth(0.2),
153 fCutTM(0.15),
154 fJetEtaMin(-.5),
155 fJetEtaMax(.5),
156 fNevents(0),
157 fTindex(0),
158 fJetPtMin(20.),
159 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
160 fJetPtFractionMin(0.5),
161 fNMatchJets(4),
162 fMatchMaxDist(0.8),
163 fKeepJets(kFALSE),
164 fkNbranches(2),
165 fkEvtClasses(12),
166 fOutputList(0x0),
167 fHistEvtSelection(0x0),
168 fh1JetEntries(0x0),
169 fh2Circularity(0x0),
170 fhnJetTM(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
184 AliAnalysisTaskJetAntenna::~AliAnalysisTaskJetAntenna()
185 {
186    delete fListJets[0];
187    delete fListJets[1];
188 }
189
190 void AliAnalysisTaskJetAntenna::SetBranchNames(const TString &branch1, const TString &branch2)
191 {
192    fJetBranchName[0] = branch1;
193    fJetBranchName[1] = branch2;
194 }
195
196 void 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
206 void 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);
229   Int_t nbinsJet[8]={10,9,7,9,36,10,2,10};
230   Double_t binlowJet[8]= {0,0, 0, 0,-0.5*TMath::Pi(),0,0,-0.5};
231   Double_t binupJet[8]= {100,0.9, 150,150,1.5*TMath::Pi(),1,200,9.5};
232   fhnJetTM = new THnSparseF("fhnJetTM", "fhnJetTM; cent;dr;pt_jet;pt_track;phi;circ;nc;pthard",8,nbinsJet,binlowJet,binupJet);
233   Double_t *xPt3=new Double_t[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(3,xPt3);
241
242   Double_t *xPt2=new Double_t[10];
243   xPt2[0] = 0.;
244   xPt2[1]=20;
245   xPt2[2]=40;
246   xPt2[3]=60;
247   xPt2[4]=80;
248   xPt2[5]=100;
249   xPt2[6]=120;
250   xPt2[7]=150; 
251   
252   fhnJetTM->SetBinEdges(2,xPt2);
253   fOutputList->Add(fhnJetTM);
254    delete [] xPt3;
255    delete [] xPt2;
256   // =========== Switch on Sumw2 for all histos ===========
257   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
258     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
259     if (h1){
260       h1->Sumw2();
261       continue;
262     }
263     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
264     if (hn){
265       hn->Sumw2();
266     }
267   }
268   TH1::AddDirectory(oldStatus);
269
270   PostData(1, fOutputList);
271 }
272
273 void AliAnalysisTaskJetAntenna::UserExec(Option_t *)
274 {
275
276
277   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
278     AliError("Jet branch name not set.");
279     return;
280   }
281
282   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
283   if (!fESD) {
284     AliError("ESD not available");
285     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
286   }
287   fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
288
289   static AliAODEvent* aod = 0;
290   // take all other information from the aod we take the tracks from
291   if(!aod){
292     if(!fESD)aod = fAODIn;
293     else aod = fAODOut;}
294
295
296
297   if(fNonStdFile.Length()!=0){
298     // case that we have an AOD extension we need can fetch the jets from the extended output
299     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
300     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
301     if(!fAODExtension){
302       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
303     }
304   }
305
306
307   // -- event selection --
308   fHistEvtSelection->Fill(1); // number of events before event selection
309
310
311   Bool_t selected=kTRUE;
312   selected = AliAnalysisHelperJetTasks::Selected();
313   if(!selected){
314     // no selection by the service task, we continue
315     PostData(1,fOutputList);
316     return;}
317
318
319
320   // physics selection: this is now redundant, all should appear as accepted after service task selection
321   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
322     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
323   std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
324   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
325     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
326     fHistEvtSelection->Fill(2);
327     PostData(1, fOutputList);
328     return;
329   }
330
331
332
333   // vertex selection
334   if(!aod){
335     if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
336     fHistEvtSelection->Fill(3);
337     PostData(1, fOutputList);
338   }
339   AliAODVertex* primVtx = aod->GetPrimaryVertex();
340
341   if(!primVtx){
342     if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
343     fHistEvtSelection->Fill(3);
344     PostData(1, fOutputList);
345     return;
346   }
347
348   Int_t nTracksPrim = primVtx->GetNContributors();
349   if ((nTracksPrim < fMinContribVtx) ||
350       (primVtx->GetZ() < fVtxZMin) ||
351       (primVtx->GetZ() > fVtxZMax) ){
352     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
353     fHistEvtSelection->Fill(3);
354     PostData(1, fOutputList);
355     return;
356   }
357
358
359
360   // centrality selection
361   AliCentrality *cent = 0x0;
362   Double_t centValue = 0.;
363   if(fIsPbPb){
364     if(fESD) {cent = fESD->GetCentrality();
365       if(cent) centValue = cent->GetCentralityPercentile("V0M");}
366     else     centValue=aod->GetHeader()->GetCentrality();
367
368     if(fDebug) printf("centrality: %f\n", centValue);
369     if (centValue < fCentMin || centValue > fCentMax){
370       fHistEvtSelection->Fill(4);
371       PostData(1, fOutputList);
372       return;
373     }}
374
375
376   fHistEvtSelection->Fill(0);
377   // accepted events
378   // -- end event selection --
379   // pt hard
380   Double_t pthard=0;
381   Double_t pthardbin=0;
382   if(fDoMatching){
383   pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
384   pthardbin = GetPtHardBin(pthard);}
385   // get background
386   AliAODJetEventBackground* externalBackground = 0;
387   if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
388     externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
389     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
390   }
391   if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
392     externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
393     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
394   }
395
396   if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
397     externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
398     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
399   }
400
401   Float_t rho = 0;
402
403
404   if(externalBackground)rho = externalBackground->GetBackground(0);
405
406   // fetch jets
407   TClonesArray *aodJets[2];
408   aodJets[0]=0;
409   if(fAODOut&&!aodJets[0]){
410     aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
411     aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));
412   }
413   if(fAODExtension && !aodJets[0]){
414     aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
415     aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
416   }
417   if(fAODIn&&!aodJets[0]){
418     aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
419     aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));
420   }
421
422
423
424   Int_t nT=0;
425   TList ParticleList;
426   if(!fDoMatching) nT = GetListOfTracks(&ParticleList);
427   if(fDoMatching)nT=GetListOfTracksExtra(&ParticleList);
428   if(nT<0){
429     PostData(1, fOutputList);
430     return;
431   }
432
433   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
434     fListJets[iJetType]->Clear();
435     if (!aodJets[iJetType]) continue;
436     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
437     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
438       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
439       if (jet) fListJets[iJetType]->Add(jet);
440     }
441   }
442
443 // jet matching 
444   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
445   static TArrayF aPtFraction(fListJets[0]->GetEntries());
446   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
447   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
448
449   if(fDoMatching){
450  
451
452   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
453   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
454                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
455                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);}
456
457   
458
459
460
461  
462  // loop over matched jets
463   Int_t ir = -1;  // index of matched reconstruced jet
464   Float_t fraction = -1.;
465    AliAODJet *jetmatched =0x0;
466    
467   
468    
469    
470   
471     for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
472
473     Double_t etabig=0;
474     Double_t ptbig=0;
475     Double_t areabig=0;
476     Double_t phibig=0.;
477     Double_t pxbig,pybig,pzbig;
478
479
480     AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
481     etabig  = jetbig->Eta();
482     phibig  = jetbig->Phi();
483     ptbig   = jetbig->Pt();
484     if(ptbig==0) continue;
485     areabig = jetbig->EffectiveAreaCharged();
486     if(!fDoMatching) ptbig=ptbig-rho*areabig;
487     if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
488
489     if(fSemigoodCorrect){
490       if((phibig>fHolePos-fHoleWidth) && (phibig<fHolePos+fHoleWidth)) continue;
491     }
492     pxbig=jetbig->Px();
493     pybig=jetbig->Py();
494     pzbig=jetbig->Pz();
495
496    if(fDoMatching){
497      Bool_t jetAccepted=kTRUE;
498     ir = aMatchIndex[i];
499     if(ir>=0) jetmatched = (AliAODJet*)(fListJets[1]->At(ir));
500     else    continue;
501     fraction  = aPtFraction[i];
502     // minimum fraction required
503     
504     if(fraction<fJetPtFractionMin) jetAccepted = kFALSE;
505     if(jetAccepted){
506     // jet acceptance + minimum pT check
507     if(jetmatched->Eta()>fJetEtaMax || jetmatched->Eta()<fJetEtaMin)  jetAccepted = kFALSE;
508       
509     }
510     if(!jetAccepted) continue;
511     pxbig=jetmatched->Px();
512     pybig=jetmatched->Py();
513     pzbig=jetmatched->Pz();
514     ptbig=jetmatched->Pt()-rho*jetmatched->EffectiveAreaCharged();}
515
516
517
518
519     TVector3  ppJ1(pxbig, pybig, pzbig);
520     TVector3  ppJ3(- pxbig * pzbig, - pybig * pzbig, pxbig * pxbig + pybig * pybig);
521     ppJ3.SetMag(1.);
522     TVector3  ppJ2(-pybig, pxbig, 0);
523     ppJ2.SetMag(1.);
524
525
526
527
528     Float_t mxx    = 0.;
529     Float_t myy    = 0.;
530     Float_t mxy    = 0.;
531     Int_t   nc     = 0;
532     Float_t sump2  = 0.;
533
534     for(int it = 0;it<ParticleList.GetEntries();++it){
535       AliVParticle *track = (AliVParticle*)ParticleList.At(it);
536       TVector3 pp(track->Px(), track->Py(), track->Pz());
537       Float_t phi = track->Phi();
538       Float_t eta = track->Eta();
539       Float_t pt  = track->Pt();
540       Float_t deta = eta - etabig;
541       Float_t dphi = RelativePhi(phi,phibig);
542       if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
543       Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
544       if (r < 0.4 && pt>fCutTM) {
545         //longitudinal and perpendicular component of the track pT in the
546         //local frame
547         TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
548         TVector3 pPerp = pp - pLong;
549         //projection onto the two perpendicular vectors defined above
550         Float_t ppjX = pPerp.Dot(ppJ2);
551         Float_t ppjY = pPerp.Dot(ppJ3);
552         Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
553         //components of the 2D symmetrical sphericity matrix
554         mxx += (ppjX * ppjX / ppjT);
555         myy += (ppjY * ppjY / ppjT);
556         mxy += (ppjX * ppjY / ppjT);
557         nc++;
558         sump2 += ppjT;}
559       // max pt
560       if(nc<2) continue;
561
562     } // 1st Track Loop
563
564
565     // Sphericity Matrix
566     const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
567     TMatrixDSym m0(2,ele);
568     // Find eigenvectors
569     TMatrixDSymEigen m(m0);
570     TVectorD eval(2);
571     TMatrixD evecm = m.GetEigenVectors();
572     eval  = m.GetEigenValues();
573     // Largest eigenvector
574     Int_t jev = 0;
575     if (eval[0] < eval[1]) jev = 1;
576     TVectorD evec0(2);
577     // Principle axis
578     evec0 = TMatrixDColumn(evecm, jev);
579     TVector2 evec(evec0[0], evec0[1]);
580     Float_t circ=0;
581     if(jev==1) circ=2*eval[0];
582     if(jev==0) circ=2*eval[1];
583     fh2Circularity->Fill(circ,ptbig);
584     fh1JetEntries->Fill(ptbig);
585     for (Int_t ip = 0; ip < ParticleList.GetEntries(); ip++) {
586       AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
587       TVector3 pp(track->Px(), track->Py(), track->Pz());
588       Float_t phi = track->Phi();
589       Float_t eta = track->Eta();
590       Float_t pt  = track->Pt();
591
592       Float_t deta = eta - etabig;
593       Float_t dphi = RelativePhi(phi,phibig);
594       if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
595
596       Float_t dRR = TMath::Sqrt(dphi * dphi + deta * deta);
597       TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
598       TVector3 pPerp = pp - pLong;
599       Float_t ppjX = pPerp.Dot(ppJ2);
600       Float_t ppjY = pPerp.Dot(ppJ3);
601       TVector2 vr(ppjX, ppjY) ;
602       //and this is the angle between the particle and the TM axis.
603       float phistr=evec.Phi()-vr.Phi();
604       if(phistr>2*TMath::Pi()) phistr -= 2*TMath::Pi();
605       if(phistr<-2*TMath::Pi()) phistr += 2*TMath::Pi();
606       if(phistr<-0.5*TMath::Pi()) phistr += 2*TMath::Pi();
607       if(phistr>1.5*TMath::Pi()) phistr -= 2*TMath::Pi();
608
609       double jetEntries[8] = {centValue,dRR,ptbig,pt,phistr,circ,static_cast<double>(nc),pthardbin};
610       fhnJetTM->Fill(jetEntries);
611
612     } // 2nd Track loop
613   }//jet loop
614
615   PostData(1, fOutputList);
616 }
617
618 void AliAnalysisTaskJetAntenna::Terminate(const Option_t *)
619 {
620    // Draw result to the screen
621    // Called once at the end of the query
622
623    if (!GetOutputData(1))
624    return;
625 }
626
627
628 Int_t  AliAnalysisTaskJetAntenna::GetListOfTracks(TList *list){
629
630   Int_t iCount = 0;
631   AliAODEvent *aod = 0;
632
633   if(!fESD)aod = fAODIn;
634   else aod = fAODOut;
635   if(!aod)return 0;
636
637
638     for(int it = 0;it < aod->GetNumberOfTracks();++it){
639     AliAODTrack *tr = aod->GetTrack(it);
640     Bool_t bGood = false;
641     if(fFilterType == 0)bGood = true;
642     else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
643     else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
644     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
645     if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
646     if(bGood==false) continue;
647     if (fApplySharedClusterCut) {
648       Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
649       if (frac > 0.4) continue;
650     }
651     if(TMath::Abs(tr->Eta())>0.9)continue;
652     if(tr->Pt()<0.15)continue;
653     list->Add(tr);
654     iCount++;
655     }
656   return iCount;
657 }
658
659
660
661
662 Int_t  AliAnalysisTaskJetAntenna::GetListOfTracksExtra(TList *list){
663
664   Int_t iCount = 0;
665   AliAODEvent *aod = 0;
666
667   if(!fESD)aod = fAODIn;
668   else aod = fAODOut;
669   if(!aod)return 0;
670  
671     for(int it = 0;it < aod->GetNumberOfTracks();++it){
672     AliAODTrack *tr = aod->GetTrack(it);
673     Bool_t bGood = false;
674     if(fFilterType == 0)bGood = true;
675     else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
676     else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
677     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
678     if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
679     if(bGood==false) continue;
680     if (fApplySharedClusterCut) {
681       Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
682       if (frac > 0.4) continue;
683     }
684     if(TMath::Abs(tr->Eta())>0.9)continue;
685     if(tr->Pt()<0.15)continue;
686     list->Add(tr);
687     iCount++;
688     } 
689
690
691
692
693
694      TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
695       if(!aodExtraTracks)return iCount;
696       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
697         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
698         if (!track) continue;
699         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
700         if(!trackAOD)continue;
701         Bool_t bGood = false;
702         if(fFilterType == 0)bGood = true;
703         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
704         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
705         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
706         if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
707          if (fApplySharedClusterCut) {
708            Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
709            if (frac > 0.4) continue;
710          }
711
712
713         if(TMath::Abs(trackAOD->Eta())>0.9) continue;
714         if(trackAOD->Pt()<0.15) continue;
715         if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
716         list->Add(trackAOD);
717         iCount++;
718       }
719     
720      return iCount;
721 }
722
723 Double_t AliAnalysisTaskJetAntenna::RelativePhi(Double_t mphi,Double_t vphi){
724
725   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
726   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
727   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
728   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
729   double dphi = mphi-vphi;
730   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
731   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
732   return dphi;//dphi in [-Pi, Pi]
733 }
734
735 Int_t AliAnalysisTaskJetAntenna::GetPhiBin(Double_t phi)
736 {
737     Int_t phibin=-1;
738     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
739     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
740     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
741     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
742     return phibin;
743 }
744
745 Int_t AliAnalysisTaskJetAntenna::GetPtHardBin(Double_t ptHard){
746
747    const Int_t nBins = 10;
748   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
749    
750   Int_t bin = -1;
751   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
752     bin++;
753   }
754    
755   return bin;
756 }