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