AliPhysics  59e0e03 (59e0e03)
AliEmcalJetTaggerTaskFast.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
27 #include <iostream>
28 #include <vector>
29 
30 #include <TH1.h>
31 #include <TH2.h>
32 #include <TH3.h>
33 #include <THnSparse.h>
34 #include <TKDTree.h>
35 
36 #include "AliAnalysisManager.h"
37 #include "AliEmcalJet.h"
38 #include "AliLog.h"
39 #include "AliJetContainer.h"
40 #include "AliParticleContainer.h"
41 
43 
47 
48 namespace PWGJE {
49 
50  namespace EMCALJetTasks{
51 
52  AliEmcalJetTaggerTaskFast::AliEmcalJetTaggerTaskFast() :
53  AliAnalysisTaskEmcalJet("AliEmcalJetTaggerTaskFast", kTRUE),
54  fJetTaggingType(kTag),
55  fJetTaggingMethod(kGeo),
56  fContainerBase(0),
57  fContainerTag(1),
58  fSpecPartContTag(-1),
59  fMinFractionShared(0),
60  fUseSumw2(0),
61  fMatchingDone(0),
62  fTypeAcc(kLimitBaseTagEtaPhi),
63  fMaxDist(0.3),
64  fInit(kFALSE),
65  fh3PtJet1VsDeltaEtaDeltaPhi(nullptr),
66  fh2PtJet1VsDeltaR(nullptr),
67  fh2PtJet2VsFraction(nullptr),
68  fh2PtJet1VsLeadPtAllSel(nullptr),
69  fh2PtJet1VsLeadPtTagged(nullptr),
70  fh2PtJet1VsPtJet2(nullptr),
71  fh2PtJet2VsRelPt(nullptr),
72  fh3PtJetDEtaDPhiConst(nullptr),
73  fh3PtJetAreaDRConst(nullptr),
74  fNAccJets(nullptr)
75 #ifdef JETTAGGERFAST_TEST
76  , fIndexErrorRateBase(nullptr)
77  , fIndexErrorRateTag(nullptr)
78  , fContainerErrorRateBase(nullptr)
79  , fContainerErrorRateTag(nullptr)
80 #endif
81  {
89 
90  for (Int_t i = 0; i < fNcentBins; i++) {
92  fh2PtJet1VsDeltaR[i] = 0;
93  fh2PtJet2VsFraction[i] = 0;
96  fh2PtJet1VsPtJet2[i] = 0;
97  fh2PtJet2VsRelPt[i] = 0;
98  }
99 
101  }
102 
104  AliAnalysisTaskEmcalJet(name, kTRUE),
107  fContainerBase(0),
108  fContainerTag(1),
109  fSpecPartContTag(-1),
111  fUseSumw2(0),
112  fMatchingDone(0),
114  fMaxDist(0.3),
115  fInit(kFALSE),
126 #ifdef JETTAGGERFAST_TEST
127  , fIndexErrorRateBase(nullptr)
128  , fIndexErrorRateTag(nullptr)
129  , fContainerErrorRateBase(nullptr)
130  , fContainerErrorRateTag(nullptr)
131 #endif
132  {
133 
141 
142  for (Int_t i = 0; i < fNcentBins; i++) {
144  fh2PtJet1VsDeltaR[i] = 0;
145  fh2PtJet2VsFraction[i] = 0;
148  fh2PtJet1VsPtJet2[i] = 0;
149  fh2PtJet2VsRelPt[i] = 0;
150  }
151 
153 
154  }
155 
158 
159  Bool_t oldStatus = TH1::AddDirectoryStatus();
160  TH1::AddDirectory(kFALSE);
161 
162  const Int_t nBinsPt = 40;
163  const Int_t nBinsDPhi = 72;
164  const Int_t nBinsDEta = 100;
165  const Int_t nBinsDR = 50;
166  const Int_t nBinsFraction = 101;
167 
168  const Double_t minPt = -50.;
169  const Double_t maxPt = 150.;
170  const Double_t minDPhi = -0.5;
171  const Double_t maxDPhi = 0.5;
172  const Double_t minDEta = -0.5;
173  const Double_t maxDEta = 0.5;
174  const Double_t minDR = 0.;
175  const Double_t maxDR = 0.5;
176  const Double_t minFraction = -0.005;
177  const Double_t maxFraction = 1.005;
178 
179  TString histName = "";
180  TString histTitle = "";
181  for (Int_t i = 0; i < fNcentBins; i++) {
182 
183  histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
184  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
185  fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
186  fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
187 
188  histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
189  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
190  fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
191  fOutput->Add(fh2PtJet1VsDeltaR[i]);
192 
193  histName = TString::Format("fh2PtJet2VsFraction_%d",i);
194  histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
195  fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
196  fOutput->Add(fh2PtJet2VsFraction[i]);
197 
198  histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
199  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
200  fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
201  fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
202 
203  histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
204  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
205  fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
206  fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
207 
208  histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
209  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
210  fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
211  fOutput->Add(fh2PtJet1VsPtJet2[i]);
212 
213  histName = TString::Format("fh2PtJet2VsRelPt_%d",i);
214  histTitle = TString::Format("%s;#it{p}_{T,jet2};(#it{p}_{T,jet2}-#it{p}_{T,jet1})/#it{p}_{T,jet1}",histName.Data());
215  fh2PtJet2VsRelPt[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,241,-2.41,2.41);
216  fOutput->Add(fh2PtJet2VsRelPt[i]);
217  }
218 
219  fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
221 
222  fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,50,0.,1.,50,0.,1.);
224 
225  fNAccJets = new TH1F("fNAccJets","fNAccJets;N/ev",11,-0.5, 9.5);
226  fOutput->Add(fNAccJets);
227 
228 #ifdef JETTAGGERFAST_TEST
229  fIndexErrorRateBase = new TH1F("indexErrorsBase", "Index errors nearest neighbor base jets", 1, 0.5, 1.5);
230  fIndexErrorRateTag = new TH1F("indexErrorsTag", "Index errors nearest neighbors tag jets", 1, 0.5, 1.5);
231  fContainerErrorRateBase = new TH1F("containerErrorsBase", "Matching errors container - kdtree base jets", 1, 0.5, 1.5);
232  fContainerErrorRateTag = new TH1F("containerErrorsTag", "Matching errors container - kdtree tag jets", 1, 0.5, 1.5);
233  fOutput->Add(fIndexErrorRateBase);
234  fOutput->Add(fIndexErrorRateTag);
235  fOutput->Add(fContainerErrorRateBase);
236  fOutput->Add(fContainerErrorRateTag);
237 #endif
238 
239  if(fUseSumw2) {
240  // =========== Switch on Sumw2 for all histos ===========
241  for(auto it : *fOutput){
242  TH1 *h1 = dynamic_cast<TH1*>(it);
243  if (h1){
244  h1->Sumw2();
245  continue;
246  }
247  THnSparse *hn = dynamic_cast<THnSparse*>(it);
248  if(hn)hn->Sumw2();
249  }
250  }
251 
252  TH1::AddDirectory(oldStatus);
253 
254  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
255  }
256 
258 
259  if(fInit) return;
260 
263  if(!cont1 || !cont2) AliError("Missing jet container");
264 
265  // when full azimuth, don't do anything
266  Double_t phiMin1 = cont1->GetJetPhiMin(), phiMin2 = cont2->GetJetPhiMin();
267  Bool_t isZeroTwoPi1 = kFALSE;
268  //check only one side of phi, since the upper bound is not well defined
269  if(phiMin1 > -1.e-6 && phiMin1 < 1.e-6) isZeroTwoPi1 = kTRUE;
270  Bool_t isZeroTwoPi2 = kFALSE;
271  if(phiMin2 > -1.e-6 && phiMin2 < 1.e-6) isZeroTwoPi2 = kTRUE;
272 
273  switch(fTypeAcc){
274  case kNoLimit: break;
275  case kLimitTagEta:
276  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
277  break;
278  case kLimitTagEtaPhi:
279  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
280  if(!isZeroTwoPi2) cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
281  break;
282  case kLimitBaseTagEtaPhi:
283  cont1->SetJetEtaLimits(cont1->GetJetEtaMin()-0.1,cont1->GetJetEtaMax()+0.1);
284  if(!isZeroTwoPi1) cont1->SetJetPhiLimits(cont1->GetJetPhiMin()-0.1,cont1->GetJetPhiMax()+0.1);
285  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
286  if(!isZeroTwoPi2) cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
287  };
288  fInit = kTRUE;
289  return;
290  }
291 
293  {
294  Init();
296  *contTag = GetJetContainer(fContainerTag);
297 
298  ResetTagging(*contBase);
299  ResetTagging(*contTag);
300 
301  fMatchingDone = MatchJetsGeo(*contBase, *contTag, fMaxDist);
302 
303  return kTRUE;
304  }
305 
307  {
308  // Fill histograms.
309 
310  AliEmcalJet *jet1 = NULL;
312  if(!jetCont) return kFALSE;
313  jetCont->ResetCurrentID();
314  Int_t count = 0;
315  while((jet1 = jetCont->GetNextAcceptJet())) {
316  count++;
317  Double_t ptJet1 = jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
318  fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
319 
320  //fill histo with angle between jet axis and constituents
321  for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
322  AliVParticle *vp = static_cast<AliVParticle*>(jet1->Track(icc));
323  if(!vp) continue;
324  Double_t dEta = jet1->Eta()-vp->Eta();
325  Double_t dPhi = jet1->Phi()-vp->Phi();
326  if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
327  if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
328  fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);
329 
330  Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
331  fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
332  }
333 
334  if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
335  continue;
336 
337  AliEmcalJet *jet2 = NULL;
338  if(fJetTaggingType==kTag) jet2 = jet1->GetTaggedJet();
339  if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
340  if(!jet2) continue;
341 
342  Double_t ptJet2 = jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
343 
344  Double_t fraction = -2;
345  if(fSpecPartContTag > -1) fraction = jetCont->GetFractionSharedPt(jet1, GetParticleContainer(fSpecPartContTag));
346  else fraction = jetCont->GetFractionSharedPt(jet1);
347 
348  fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);
349  AliDebug(5, Form("Fraction = %f, minimum = %f", fraction, fMinFractionShared));
350  //if(fJetTaggingType==kClosest) Printf("Fraction = %f, minimum = %f", fraction, fMinFractionShared);
352  continue;
353  fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
354  fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
355  if(ptJet2>0.) fh2PtJet2VsRelPt[fCentBin]->Fill(ptJet2,(ptJet1-ptJet2)/ptJet2);
356 
357  Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
358  if(dPhi>TMath::Pi())
359  dPhi -= TMath::TwoPi();
360  if(dPhi<(-1.*TMath::Pi()))
361  dPhi += TMath::TwoPi();
362 
363  fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
364  fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,jet1->DeltaR(jet2));
365  }
366  fNAccJets->Fill(count);
367  return kTRUE;
368  }
369 
371  for(auto j : c.all()){
372  switch(fJetTaggingType){
373  case kClosest:
374  j->ResetMatching();
375  break;
376  case kTag:
377  j->SetTaggedJet(nullptr);
378  j->SetTagStatus(-1);
379  };
380  }
381  }
382 
384  const Int_t kNacceptedBase = contBase.GetNAcceptedJets(),
385  kNacceptedTag = contTag.GetNAcceptedJets();
386  if(!(kNacceptedBase && kNacceptedTag)) return false;
387 
388  // Build kd-trees
389  TArrayD etaBase(kNacceptedBase), phiBase(kNacceptedBase),
390  etaTag(kNacceptedTag), phiTag(kNacceptedTag);
391  std::vector<AliEmcalJet *> jetsBase(kNacceptedBase), jetsTag(kNacceptedTag); // the storages are needed later for applying the tagging, in order to avoid multiple occurrence of jet selection
392 
393  int countBase(0), countTag(0);
394  for(auto jb : contBase.accepted()) {
395  etaBase[countBase] = jb->Eta();
396  phiBase[countBase] = jb->Phi();
397  jetsBase[countBase] = jb;
398  countBase++;
399  }
400  for(auto jt : contTag.accepted()) {
401  etaTag[countTag] = jt->Eta();
402  phiTag[countTag] = jt->Phi();
403  jetsTag[countTag] = jt;
404  countTag++;
405  }
406  TKDTreeID treeBase(etaBase.GetSize(), 2, 1), treeTag(etaTag.GetSize(), 2, 1);
407  treeBase.SetData(0, etaBase.GetArray());
408  treeBase.SetData(1, phiBase.GetArray());
409  treeBase.Build();
410  treeTag.SetData(0, etaTag.GetArray());
411  treeTag.SetData(1, phiTag.GetArray());
412  treeTag.Build();
413 
414  TArrayI faMatchIndexTag(kNacceptedBase), faMatchIndexBase(kNacceptedTag);
415  faMatchIndexBase.Reset(-1);
416  faMatchIndexTag.Reset(-1);
417 
418  // find the closest distance to the full jet
419  countBase = 0;
420  for(auto j : contBase.accepted()) {
421  Double_t point[2] = {j->Eta(), j->Phi()};
422  Int_t index(-1); Double_t distance(-1);
423  treeTag.FindNearestNeighbors(point, 1, &index, &distance);
424  // test whether indices are matching:
425  if(index >= 0 && distance < maxDist){
426  AliDebugStream(1) << "Found closest tag jet for " << countBase << " with match index " << index << " and distance " << distance << std::endl;
427  faMatchIndexTag[countBase]=index;
428  } else {
429  AliDebugStream(1) << "Not found closest tag jet for " << countBase << ", distance to closest " << distance << std::endl;
430  }
431 
432 #ifdef JETTAGGERFAST_TEST
433  if(index>-1){
434  Double_t distanceTest(-1);
435  distanceTest = TMath::Sqrt(TMath::Power(etaTag[index] - j->Eta(), 2) + TMath::Power(phiTag[index] - j->Phi(), 2));
436  if(TMath::Abs(distanceTest - distance) > DBL_EPSILON){
437  AliDebugStream(1) << "Mismatch in distance from tag jet with index from tree: " << distanceTest << ", distance from tree " << distance << std::endl;
438  fIndexErrorRateBase->Fill(1);
439  }
440  }
441 #endif
442 
443  countBase++;
444  }
445 
446  // other way around
447  countTag = 0;
448  for(auto j : contTag.accepted()){
449  Double_t point[2] = {j->Eta(), j->Phi()};
450  Int_t index(-1); Double_t distance(-1);
451  treeBase.FindNearestNeighbors(point, 1, &index, &distance);
452  if(index >= 0 && distance < maxDist){
453  AliDebugStream(1) << "Found closest base jet for " << countBase << " with match index " << index << " and distance " << distance << std::endl;
454  faMatchIndexBase[countTag]=index;
455  } else {
456  AliDebugStream(1) << "Not found closest tag jet for " << countBase << ", distance to closest " << distance << std::endl;
457  }
458 
459 #ifdef JETTAGGERFAST_TEST
460  if(index>-1){
461  Double_t distanceTest(-1);
462  distanceTest = TMath::Sqrt(TMath::Power(etaBase[index] - j->Eta(), 2) + TMath::Power(phiBase[index] - j->Phi(), 2));
463  if(TMath::Abs(distanceTest - distance) > DBL_EPSILON){
464  AliDebugStream(1) << "Mismatch in distance from base jet with index from tree: " << distanceTest << ", distance from tree " << distance << std::endl;
465  fIndexErrorRateTag->Fill(1);
466  }
467  }
468 #endif
469 
470  countTag++;
471  }
472 
473  // check for "true" correlations
474  // these are pairs where the base jet is the closest to the tag jet and vice versa
475  // As the lists are linear a loop over the outer base jet is sufficient.
476  AliDebugStream(1) << "Starting true jet loop: nbase(" << kNacceptedBase << "), ntag(" << kNacceptedTag << ")\n";
477  for(int ibase = 0; ibase < kNacceptedBase; ibase++) {
478  AliDebugStream(2) << "base jet " << ibase << ": match index in tag jet container " << faMatchIndexTag[ibase] << "\n";
479  if(faMatchIndexTag[ibase] > -1){
480  AliDebugStream(2) << "tag jet " << faMatchIndexTag[ibase] << ": matched base jet " << faMatchIndexBase[faMatchIndexTag[ibase]] << "\n";
481  }
482  if(faMatchIndexTag[ibase] > -1 && faMatchIndexBase[faMatchIndexTag[ibase]] == ibase) {
483  AliDebugStream(2) << "found a true match \n";
484  AliEmcalJet *jetBase = jetsBase[ibase],
485  *jetTag = jetsTag[faMatchIndexTag[ibase]];
486  if(jetBase && jetTag) {
487 #ifdef JETTAGGERFAST_TEST
488  if(TMath::Abs(etaBase[ibase] - jetBase->Eta()) > DBL_EPSILON || TMath::Abs(phiBase[ibase] - jetBase->Phi()) > DBL_EPSILON){
489  AliErrorStream() << "Selected incorrect base jet for tagging : eta test(" << jetBase->Eta() << ")/true(" << etaBase[ibase]
490  << "), phi test(" << jetBase->Phi() << ")/true(" << phiBase[ibase] << ")\n";
491  fContainerErrorRateBase->Fill(1);
492  }
493  if(TMath::Abs(etaTag[faMatchIndexTag[ibase]] - jetTag->Eta()) > DBL_EPSILON || TMath::Abs(phiTag[faMatchIndexTag[ibase]] - jetTag->Phi()) > DBL_EPSILON){
494  AliErrorStream() << "Selected incorrect tag jet for tagging : eta test(" << jetTag->Eta() << ")/true(" << etaTag[faMatchIndexTag[ibase]]
495  << "), phi test(" << jetTag->Phi() << ")/true(" << phiTag[faMatchIndexTag[ibase]] << ")\n";
496  fContainerErrorRateTag->Fill(1);
497  }
498 #endif
499  // Test if the position of the jets correp
500  Double_t dR = jetBase->DeltaR(jetTag);
501  switch(fJetTaggingType){
502  case kTag:
503  jetBase->SetTaggedJet(jetTag);
504  jetBase->SetTagStatus(1);
505 
506  jetTag->SetTaggedJet(jetBase);
507  jetTag->SetTagStatus(1);
508  break;
509  case kClosest:
510  jetBase->SetClosestJet(jetTag,dR);
511  jetTag->SetClosestJet(jetBase,dR);
512  break;
513  };
514  }
515  }
516  }
517  return kTRUE;
518  }
519 
521  return GetDeltaPhi(jet1->Phi(),jet2->Phi());
522  }
523 
525  Double_t dPhi = phi1-phi2;
526  if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
527  if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
528 
529  return dPhi;
530  }
531 
533  const char * njetsTag,
534  const Double_t R,
535  const char * nrhoBase,
536  const char * nrhoTag,
537  const char * ntracks,
538  const char * nclusters,
539  const char * type,
540  const char * CentEst,
541  Int_t pSel,
542  const char * trigClass){
543 
544  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
545  if (!mgr) {
546  std::cerr << "E-AddEmcalJetTaggerTaskFast: No analysis manager found.\n";
547  return nullptr;
548  }
549  // Check the analysis type using the event handlers connected to the analysis manager.
550  //==============================================================================
551  if (!mgr->GetInputEventHandler())
552  {
553  std::cerr << "E-AddEmcalJetTaggerTaskFast: This task requires an input event handler\n";
554  return NULL;
555  }
556 
557  TString wagonName = Form("JetTagger_%s_%s_TC%s",njetsBase,njetsTag,trigClass);
558 
559  //Configure jet tagger task
561 
562  task->SetNCentBins(4);
563  //task->SetVzRange(-10.,10.);
564 
565  AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
566  AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
567 
568  task->SetJetContainerBase(0);
569  task->SetJetContainerTag(1);
570 
571  TString strType(type);
572  AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
573  if(jetContBase) {
574  jetContBase->SetRhoName(nrhoBase);
575  jetContBase->ConnectParticleContainer(trackCont);
576  jetContBase->ConnectClusterContainer(clusterCont);
577  jetContBase->SetMaxTrackPt(10000.);
578  }
579 
580  AliJetContainer *jetContTag = task->AddJetContainer(njetsTag,"TPC",R);
581  if(jetContTag) {
582  jetContTag->SetRhoName(nrhoTag);
583  jetContTag->ConnectParticleContainer(trackCont);
584  jetContTag->ConnectClusterContainer(clusterCont);
585  jetContTag->SetMaxTrackPt(10000.);
586  }
587  for(Int_t i=0; i<2; i++) {
588  task->SetPercAreaCut(0.6, i); //keep?
589  }
590  task->SetCentralityEstimator(CentEst);
591  task->SelectCollisionCandidates(pSel);
592  task->SetUseAliAnaUtils(kFALSE);
593 
594  mgr->AddTask(task);
595 
596  //Connnect input
597  mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
598 
599  //Connect output
600  TString contName(wagonName);
601  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
602  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
603  mgr->ConnectOutput(task,1,coutput1);
604 
605  return task;
606  }
607 
608 
609  }
610 }
611 
void SetTagStatus(Int_t i)
Definition: AliEmcalJet.h:336
void SetPercAreaCut(Float_t p, Int_t c=0)
Double_t Area() const
Definition: AliEmcalJet.h:130
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
AliEmcalJet * GetTaggedJet() const
Definition: AliEmcalJet.h:337
Definition: External.C:260
Definition: External.C:236
Double_t GetJetEtaMin() const
void SetTaggedJet(AliEmcalJet *j)
Definition: AliEmcalJet.h:335
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
Int_t GetTagStatus() const
Definition: AliEmcalJet.h:338
Definition: External.C:244
TH2 ** fh2PtJet1VsPtJet2
! pT of base jet vs tagged jet
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t GetJetPhiMax() const
Double_t Phi() const
Definition: AliEmcalJet.h:117
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
Adding 0.1 in eta and phi limits for both base and tag jets.
Double_t GetJetEtaMax() const
AliVParticle * Track(Int_t idx) const
Double_t fMaxDist
distance allowed for two jets to match
TCanvas * c
Definition: TestFitELoss.C:172
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Int_t fCentBin
!event centrality bin
Double_t GetDeltaPhi(const AliEmcalJet *jet1, const AliEmcalJet *jet2)
Calculate azimuthal angle between the axes of the jets.
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
void SetRhoName(const char *n)
Bool_t fInit
true when the containers are initialized
Double_t GetJetPhiMin() const
Bool_t FillHistograms()
Filling QA histograms monitoring the quality of the matching.
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
static AliEmcalJetTaggerTaskFast * AddTaskJetTaggerFast(const char *njetsBase, const char *njetsTag, const Double_t R, const char *nrhoBase, const char *nrhoTag, const char *ntracks, const char *nclusters, const char *type, const char *CentEst, Int_t pSel, const char *trigClass)
Factory creating new jet matching task.
void SetJetPhiLimits(Float_t min, Float_t max)
TH3 * fh3PtJetDEtaDPhiConst
! jet vs delta eta vs delta phi of constituents
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
Int_t fNcentBins
how many centrality bins
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:155
void Init()
Initializing the matching task by defining the accepted range on both jet containers.
Bool_t fMatchingDone
flag to indicate if matching is done or not
AliEmcalJet * GetNextAcceptJet()
Double_t DeltaR(const AliVParticle *part) const
virtual void SetNCentBins(Int_t n)
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
TH2 ** fh2PtJet1VsLeadPtAllSel
! all jets after std selection
Double_t GetRhoVal(Int_t i=0) const
Bool_t fUseSumw2
activate sumw2 for output histograms
TH2 ** fh2PtJet2VsRelPt
! pT of tagged jet vs pt base jet / pt tagged jet
AliEmcalList * fOutput
!output list
void SetClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:323
void ResetTagging(const AliJetContainer &cont) const
Reset tagging for all jets in jet container.
Definition: External.C:220
JetTaggingMethod fJetTaggingMethod
jet matching method
TH3 * fh3PtJetAreaDRConst
! jet vs Area vs delta R of constituents
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
AcceptanceType fTypeAcc
acceptance cut for the jet containers, see method MatchJetsGeo in .cxx for possibilities ...
void UserCreateOutputObjects()
Main initialization function on the worker.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Double_t fMinFractionShared
only fill histos for jets if shared fraction larger than X
void SetCentralityEstimator(const char *c)
Int_t fSpecPartContTag
particle container optionally used in AliJetContainer::GetFractionSharedPt(). Set only if needed...
void ConnectClusterContainer(AliClusterContainer *c)
void SetMaxTrackPt(Float_t b)
void SetJetEtaLimits(Float_t min, Float_t max)
Container structure for EMCAL clusters.
No Additional limit compared to jet containers.
Container for jet within the EMCAL jet framework.
Bool_t MatchJetsGeo(AliJetContainer &contBase, AliJetContainer &contTag, Float_t maxDist=0.3) const
Match the full jets to the corresponding charged jets.
Definition: External.C:196
Fast jet tagger for geometric matching of jets.
const AliJetIterableContainer all() const