2 // Jet tagger analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.h"
24 #include "AliEmcalParticle.h"
25 #include "AliMCEvent.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliAODMCHeader.h"
28 #include "AliMCEvent.h"
29 #include "AliAnalysisManager.h"
30 #include "AliJetContainer.h"
31 #include "AliParticleContainer.h"
33 #include "AliAODEvent.h"
34 #include "AliESDEvent.h"
36 #include "AliAnalysisTaskEmcalJetTagger.h"
38 ClassImp(AliAnalysisTaskEmcalJetTagger)
40 //________________________________________________________________________
41 AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
42 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTagger", kTRUE),
43 fJetTaggingType(kTag),
44 fJetTaggingMethod(kGeo),
47 fMinFractionShared(0),
49 fh3PtJet1VsDeltaEtaDeltaPhi(0),
51 fh2PtJet2VsFraction(0),
52 fh2PtJet1VsLeadPtAllSel(0),
53 fh2PtJet1VsLeadPtTagged(0),
56 fh3PtJetDEtaDPhiConst(0),
58 fh3PtJetAreaDRConst(0)
60 // Default constructor.
62 fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
63 fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
64 fh2PtJet2VsFraction = new TH2F*[fNcentBins];
65 fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
66 fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
67 fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
68 fh2PtJet2VsRelPt = new TH2F*[fNcentBins];
70 for (Int_t i = 0; i < fNcentBins; i++) {
71 fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
72 fh2PtJet1VsDeltaR[i] = 0;
73 fh2PtJet2VsFraction[i] = 0;
74 fh2PtJet1VsLeadPtAllSel[i] = 0;
75 fh2PtJet1VsLeadPtTagged[i] = 0;
76 fh2PtJet1VsPtJet2[i] = 0;
77 fh2PtJet2VsRelPt[i] = 0;
80 SetMakeGeneralHistograms(kTRUE);
83 //________________________________________________________________________
84 AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
85 AliAnalysisTaskEmcalJet(name, kTRUE),
86 fJetTaggingType(kTag),
87 fJetTaggingMethod(kGeo),
90 fMinFractionShared(0),
92 fh3PtJet1VsDeltaEtaDeltaPhi(0),
94 fh2PtJet2VsFraction(0),
95 fh2PtJet1VsLeadPtAllSel(0),
96 fh2PtJet1VsLeadPtTagged(0),
99 fh3PtJetDEtaDPhiConst(0),
101 fh3PtJetAreaDRConst(0)
103 // Standard constructor.
105 fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
106 fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
107 fh2PtJet2VsFraction = new TH2F*[fNcentBins];
108 fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
109 fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
110 fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
111 fh2PtJet2VsRelPt = new TH2F*[fNcentBins];
113 for (Int_t i = 0; i < fNcentBins; i++) {
114 fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
115 fh2PtJet1VsDeltaR[i] = 0;
116 fh2PtJet2VsFraction[i] = 0;
117 fh2PtJet1VsLeadPtAllSel[i] = 0;
118 fh2PtJet1VsLeadPtTagged[i] = 0;
119 fh2PtJet1VsPtJet2[i] = 0;
120 fh2PtJet2VsRelPt[i] = 0;
123 SetMakeGeneralHistograms(kTRUE);
126 //________________________________________________________________________
127 AliAnalysisTaskEmcalJetTagger::~AliAnalysisTaskEmcalJetTagger()
132 //________________________________________________________________________
133 void AliAnalysisTaskEmcalJetTagger::UserCreateOutputObjects()
135 // Create user output.
137 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
139 Bool_t oldStatus = TH1::AddDirectoryStatus();
140 TH1::AddDirectory(kFALSE);
142 const Int_t nBinsPt = 250;
143 const Int_t nBinsDPhi = 100;
144 const Int_t nBinsDEta = 100;
145 const Int_t nBinsDR = 50;
146 const Int_t nBinsFraction = 101;
148 const Double_t minPt = -50.;
149 const Double_t maxPt = 200.;
150 const Double_t minDPhi = -0.5;
151 const Double_t maxDPhi = 0.5;
152 const Double_t minDEta = -0.5;
153 const Double_t maxDEta = 0.5;
154 const Double_t minDR = 0.;
155 const Double_t maxDR = 0.5;
156 const Double_t minFraction = -0.005;
157 const Double_t maxFraction = 1.005;
159 TString histName = "";
160 TString histTitle = "";
161 for (Int_t i = 0; i < fNcentBins; i++) {
163 histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
164 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
165 fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
166 fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
168 histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
169 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
170 fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
171 fOutput->Add(fh2PtJet1VsDeltaR[i]);
173 histName = TString::Format("fh2PtJet2VsFraction_%d",i);
174 histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
175 fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
176 fOutput->Add(fh2PtJet2VsFraction[i]);
178 histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
179 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
180 fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
181 fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
183 histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
184 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
185 fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
186 fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
188 histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
189 histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
190 fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
191 fOutput->Add(fh2PtJet1VsPtJet2[i]);
193 histName = TString::Format("fh2PtJet2VsRelPt_%d",i);
194 histTitle = TString::Format("%s;#it{p}_{T,jet2};(#it{p}_{T,jet2}-#it{p}_{T,jet1})/#it{p}_{T,jet1}",histName.Data());
195 fh2PtJet2VsRelPt[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,241,-2.41,2.41);
196 fOutput->Add(fh2PtJet2VsRelPt[i]);
199 fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
200 fOutput->Add(fh3PtJetDEtaDPhiConst);
202 fh2PtJetDRConst = new TH2F("fh2PtJetDRConst","fh2PtJetDRConst;pT;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.);
203 fOutput->Add(fh2PtJetDRConst);
205 fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.,100,0.,1.);
206 fOutput->Add(fh3PtJetAreaDRConst);
208 // =========== Switch on Sumw2 for all histos ===========
209 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
210 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
215 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
219 TH1::AddDirectory(oldStatus);
221 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
224 //________________________________________________________________________
225 Bool_t AliAnalysisTaskEmcalJetTagger::Run()
227 // Run analysis code here, if needed. It will be executed before FillHistograms().
229 MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,2);
231 // if(fJetTaggingMethod==kFraction)
236 //________________________________________________________________________
237 Bool_t AliAnalysisTaskEmcalJetTagger::FillHistograms()
241 AliEmcalJet *jet1 = NULL;
242 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
243 if(!jetCont) return kFALSE;
244 jetCont->ResetCurrentID();
245 while((jet1 = jetCont->GetNextAcceptJet())) {
246 Double_t ptJet1 = jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
247 fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
249 //fill histo with angle between jet axis and constituents
250 for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
251 AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(icc, jetCont->GetParticleContainer()->GetArray()));//fTracks));
253 Double_t dEta = jet1->Eta()-vp->Eta();
254 Double_t dPhi = jet1->Phi()-vp->Phi();
255 if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
256 if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
257 fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);
259 Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
260 fh2PtJetDRConst->Fill(ptJet1,dR);
261 fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
264 if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
267 AliEmcalJet *jet2 = NULL;
268 if(fJetTaggingType==kTag) jet2 = jet1->GetTaggedJet();
269 if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
272 Double_t ptJet2 = jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
273 Double_t fraction = jetCont->GetFractionSharedPt(jet1);
274 fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);
276 if(fraction<fMinFractionShared && fJetTaggingType==kClosest)
279 fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
280 fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
281 if(ptJet2>0.) fh2PtJet2VsRelPt[fCentBin]->Fill(ptJet2,(ptJet1-ptJet2)/ptJet2);
283 Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
285 dPhi -= TMath::TwoPi();
286 if(dPhi<(-1.*TMath::Pi()))
287 dPhi += TMath::TwoPi();
289 fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
290 fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,GetDeltaR(jet1,jet2));
295 //________________________________________________________________________
296 void AliAnalysisTaskEmcalJetTagger::ResetTagging(const Int_t c) {
298 //Reset tagging of container c
300 for(int i = 0;i<GetNJets(c);i++){
301 AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetJetFromArray(i, c));
302 if(fJetTaggingType==kClosest)
303 jet->ResetMatching();
304 else if(fJetTaggingType==kTag) {
305 jet->SetTaggedJet(0x0);
306 jet->SetTagStatus(-1);
311 //________________________________________________________________________
312 void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
313 Int_t iDebug, Float_t maxDist, Int_t type, Bool_t bReset) {
316 // Match the full jets to the corresponding charged jets
317 // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
318 // type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for c2 in eta 2 = allow 0.1 more in eta and phi for c2
320 if(c1<0) c1 = fContainerBase;
321 if(c2<0) c2 = fContainerTag;
323 const Int_t nJets1 = GetNJets(c1);
324 const Int_t nJets2 = GetNJets(c2);
326 if(nJets1==0 || nJets2==0) return;
332 fMatchingDone = kFALSE;
334 TArrayI faMatchIndex1;
335 faMatchIndex1.Set(nJets2+1);
336 faMatchIndex1.Reset(-1);
338 TArrayI faMatchIndex2;
339 faMatchIndex2.Set(nJets1+1);
340 faMatchIndex2.Reset(-1);
342 static TArrayS iFlag(nJets1*nJets2);
343 if(iFlag.GetSize()<(nJets1*nJets2)){
344 iFlag.Set(nJets1*nJets2+1);
349 AliJetContainer *cont2 = GetJetContainer(c2);
351 // find the closest distance to the full jet
352 for(int i = 0;i<nJets1;i++){
354 AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
357 Float_t dist = maxDist;
359 for(int j = 0;j <nJets2; j++){
360 AliEmcalJet *jet2 = 0x0;
362 jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
364 jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
367 if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
370 if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
378 Double_t dR = GetDeltaR(jet1,jet2);
379 if(dR<dist && dR<maxDist){
384 if(faMatchIndex2[i]>=0) {
385 iFlag[i*nJets2+faMatchIndex2[i]]+=1;//j closest to i
386 if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",i,faMatchIndex2[i],dist,i*nJets2+faMatchIndex2[i],iFlag[i*nJets2+faMatchIndex2[i]]);
393 for(int j = 0;j<nJets2;j++){
394 AliEmcalJet *jet2 = 0x0;
396 jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
398 jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
401 if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
404 if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
412 Float_t dist = maxDist;
413 for(int i = 0;i<nJets1;i++){
414 AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
418 Double_t dR = GetDeltaR(jet1,jet2);
419 if(dR<dist && dR<maxDist){
424 if(faMatchIndex1[j]>=0) {
425 iFlag[faMatchIndex1[j]*nJets2+j]+=2;//i closest to j
426 if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faMatchIndex1[j],j,dist,faMatchIndex1[j]*nJets2+j,iFlag[faMatchIndex1[j]*nJets2+j]);
431 // check for "true" correlations
432 for(int i = 0;i<nJets1;i++){
433 AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetJetFromArray(i, c1));
434 for(int j = 0;j<nJets2;j++){
435 AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
436 AliDebug(11,Form("%s: Flag[%d][%d] %d ",GetName(),i,j,iFlag[i*nJets2+j]));
438 // we have a uniqe correlation
439 if(iFlag[i*nJets2+j]==3){
441 Double_t dR = GetDeltaR(jet1,jet2);
442 if(iDebug>1) Printf("closest jets %d %d dR = %f",j,i,dR);
444 if(fJetTaggingType==kTag) {
445 jet1->SetTaggedJet(jet2);
446 jet1->SetTagStatus(1);
448 jet2->SetTaggedJet(jet1);
449 jet2->SetTagStatus(1);
451 else if(fJetTaggingType==kClosest) {
452 jet1->SetClosestJet(jet2,dR);
453 jet2->SetClosestJet(jet1,dR);
458 fMatchingDone = kTRUE;
461 //________________________________________________________________________
462 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
464 // Helper function to calculate the distance between two jets
467 Double_t dPhi = jet1->Phi() - jet2->Phi();
469 dPhi -= TMath::TwoPi();
470 if(dPhi<(-1.*TMath::Pi()))
471 dPhi += TMath::TwoPi();
472 Double_t dEta = jet1->Eta() - jet2->Eta();
473 Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
477 //________________________________________________________________________
478 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
480 // Calculate azimuthal angle between the axises of the jets
483 return GetDeltaPhi(jet1->Phi(),jet2->Phi());
487 //________________________________________________________________________
488 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(Double_t phi1,Double_t phi2) {
490 // Calculate azimuthal angle between the axises of the jets
493 Double_t dPhi = phi1-phi2;
495 if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
496 if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
501 //________________________________________________________________________
502 Double_t AliAnalysisTaskEmcalJetTagger::GetFractionSharedPt(const AliEmcalJet *jet1, const AliEmcalJet *jet2) const
505 // Get fraction of shared pT between matched full and charged jet
506 // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
509 Double_t fraction = 0.;
510 Double_t jetPt2 = jet2->Pt();
514 AliDebug(2,Form("%s: nConstituents_jet1: %d, nConstituents_jet2: %d, tracks_jet1: %d tracks_jet2: %d clusters_jet1: %d clusters_jet2: %d",GetName(),jet1->GetNumberOfConstituents(),jet2->GetNumberOfConstituents(),jet1->GetNumberOfTracks(),jet2->GetNumberOfTracks(),jet1->GetNumberOfClusters(),jet2->GetNumberOfClusters()));
517 AliVParticle *vpf = 0x0;
519 for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
520 Int_t idx = (Int_t)jet2->TrackAt(icc);
522 for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
523 if(idx == jet1->TrackAt(icf) && iFound==0 ) {
525 vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fTracks));
526 if(vpf) sumPt += vpf->Pt();
531 fraction = sumPt/jetPt2;
532 } else fraction = -1.;
536 //________________________________________________________________________
537 Bool_t AliAnalysisTaskEmcalJetTagger::RetrieveEventObjects() {
539 // retrieve event objects
542 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
549 //_______________________________________________________________________
550 void AliAnalysisTaskEmcalJetTagger::Terminate(Option_t *)
552 // Called once at the end of the analysis.