AliPhysics  master (3d17d9d)
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  {
83  }
84 
86  AliAnalysisTaskEmcalJet(name, kTRUE),
89  fContainerBase(0),
90  fContainerTag(1),
91  fSpecPartContTag(-1),
93  fUseSumw2(0),
94  fMatchingDone(0),
96  fMaxDist(0.3),
97  fInit(kFALSE),
108 #ifdef JETTAGGERFAST_TEST
109  , fIndexErrorRateBase(nullptr)
110  , fIndexErrorRateTag(nullptr)
111  , fContainerErrorRateBase(nullptr)
112  , fContainerErrorRateTag(nullptr)
113 #endif
114  {
116  }
117 
120 
121  // Notify the user to be careful.
122  AliErrorStream() << "This task isn't yet validated. Please use the standard AliAnalysisTaskEmcalJetTagger.\n";
123 
124  Bool_t oldStatus = TH1::AddDirectoryStatus();
125  TH1::AddDirectory(kFALSE);
126 
127  const Int_t nBinsPt = 40;
128  const Int_t nBinsDPhi = 72;
129  const Int_t nBinsDEta = 100;
130  const Int_t nBinsDR = 50;
131  const Int_t nBinsFraction = 101;
132 
133  const Double_t minPt = -50.;
134  const Double_t maxPt = 150.;
135  const Double_t minDPhi = -0.5;
136  const Double_t maxDPhi = 0.5;
137  const Double_t minDEta = -0.5;
138  const Double_t maxDEta = 0.5;
139  const Double_t minDR = 0.;
140  const Double_t maxDR = 0.5;
141  const Double_t minFraction = -0.005;
142  const Double_t maxFraction = 1.005;
143 
144  // Prepare histograms
152 
153  for (Int_t i = 0; i < fNcentBins; i++) {
155  fh2PtJet1VsDeltaR[i] = 0;
156  fh2PtJet2VsFraction[i] = 0;
159  fh2PtJet1VsPtJet2[i] = 0;
160  fh2PtJet2VsRelPt[i] = 0;
161  }
162 
163  TString histName = "";
164  TString histTitle = "";
165  for (Int_t i = 0; i < fNcentBins; i++) {
166 
167  histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
168  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
169  fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
170  fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
171 
172  histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
173  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
174  fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
175  fOutput->Add(fh2PtJet1VsDeltaR[i]);
176 
177  histName = TString::Format("fh2PtJet2VsFraction_%d",i);
178  histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
179  fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
180  fOutput->Add(fh2PtJet2VsFraction[i]);
181 
182  histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
183  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
184  fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
185  fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
186 
187  histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
188  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
189  fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
190  fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
191 
192  histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
193  histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
194  fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
195  fOutput->Add(fh2PtJet1VsPtJet2[i]);
196 
197  histName = TString::Format("fh2PtJet2VsRelPt_%d",i);
198  histTitle = TString::Format("%s;#it{p}_{T,jet2};(#it{p}_{T,jet2}-#it{p}_{T,jet1})/#it{p}_{T,jet1}",histName.Data());
199  fh2PtJet2VsRelPt[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,241,-2.41,2.41);
200  fOutput->Add(fh2PtJet2VsRelPt[i]);
201  }
202 
203  fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
205 
206  fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,50,0.,1.,50,0.,1.);
208 
209  fNAccJets = new TH1F("fNAccJets","fNAccJets;N/ev",10,-0.5, 9.5);
210  fOutput->Add(fNAccJets);
211 
212 #ifdef JETTAGGERFAST_TEST
213  fIndexErrorRateBase = new TH1F("indexErrorsBase", "Index errors nearest neighbor base jets", 1, 0.5, 1.5);
214  fIndexErrorRateTag = new TH1F("indexErrorsTag", "Index errors nearest neighbors tag jets", 1, 0.5, 1.5);
215  fContainerErrorRateBase = new TH1F("containerErrorsBase", "Matching errors container - kdtree base jets", 1, 0.5, 1.5);
216  fContainerErrorRateTag = new TH1F("containerErrorsTag", "Matching errors container - kdtree tag jets", 1, 0.5, 1.5);
217  fOutput->Add(fIndexErrorRateBase);
218  fOutput->Add(fIndexErrorRateTag);
219  fOutput->Add(fContainerErrorRateBase);
220  fOutput->Add(fContainerErrorRateTag);
221 #endif
222 
223  if(fUseSumw2) {
224  // =========== Switch on Sumw2 for all histos ===========
225  for(auto it : *fOutput){
226  TH1 *h1 = dynamic_cast<TH1*>(it);
227  if (h1){
228  h1->Sumw2();
229  continue;
230  }
231  THnSparse *hn = dynamic_cast<THnSparse*>(it);
232  if(hn)hn->Sumw2();
233  }
234  }
235 
236  TH1::AddDirectory(oldStatus);
237 
238  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
239  }
240 
242 
243  if(fInit) return;
244 
247  if(!cont1 || !cont2) AliError("Missing jet container");
248 
249  // when full azimuth, don't do anything
250  Double_t phiMin1 = cont1->GetJetPhiMin(), phiMin2 = cont2->GetJetPhiMin();
251  Bool_t isZeroTwoPi1 = kFALSE;
252  //check only one side of phi, since the upper bound is not well defined
253  if(phiMin1 > -1.e-6 && phiMin1 < 1.e-6) isZeroTwoPi1 = kTRUE;
254  Bool_t isZeroTwoPi2 = kFALSE;
255  if(phiMin2 > -1.e-6 && phiMin2 < 1.e-6) isZeroTwoPi2 = kTRUE;
256 
257  switch(fTypeAcc){
258  case kNoLimit: break;
259  case kLimitTagEta:
260  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
261  break;
262  case kLimitTagEtaPhi:
263  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
264  if(!isZeroTwoPi2) cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
265  break;
266  case kLimitBaseTagEtaPhi:
267  cont1->SetJetEtaLimits(cont1->GetJetEtaMin()-0.1,cont1->GetJetEtaMax()+0.1);
268  if(!isZeroTwoPi1) cont1->SetJetPhiLimits(cont1->GetJetPhiMin()-0.1,cont1->GetJetPhiMax()+0.1);
269  cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
270  if(!isZeroTwoPi2) cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
271  };
272  fInit = kTRUE;
273  return;
274  }
275 
277  {
278  Init();
280  *contTag = GetJetContainer(fContainerTag);
281 
282  ResetTagging(*contBase);
283  ResetTagging(*contTag);
284 
285  fMatchingDone = MatchJetsGeo(*contBase, *contTag, fMaxDist);
286 
287  return kTRUE;
288  }
289 
291  {
292  // Fill histograms.
293 
294  AliEmcalJet *jet1 = NULL;
296  if(!jetCont) return kFALSE;
297  jetCont->ResetCurrentID();
298  Int_t count = 0;
299  while((jet1 = jetCont->GetNextAcceptJet())) {
300  count++;
301  Double_t ptJet1 = jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
302  fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
303 
304  //fill histo with angle between jet axis and constituents
305  for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
306  AliVParticle *vp = static_cast<AliVParticle*>(jet1->Track(icc));
307  if(!vp) continue;
308  Double_t dEta = jet1->Eta()-vp->Eta();
309  Double_t dPhi = jet1->Phi()-vp->Phi();
310  if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
311  if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
312  fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);
313 
314  Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
315  fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
316  }
317 
318  if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
319  continue;
320 
321  AliEmcalJet *jet2 = NULL;
322  if(fJetTaggingType==kTag) jet2 = jet1->GetTaggedJet();
323  if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
324  if(!jet2) continue;
325 
326  Double_t ptJet2 = jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
327 
328  Double_t fraction = -2;
329  if(fSpecPartContTag > -1) fraction = jetCont->GetFractionSharedPt(jet1, GetParticleContainer(fSpecPartContTag));
330  else fraction = jetCont->GetFractionSharedPt(jet1);
331 
332  fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);
333  AliDebug(5, Form("Fraction = %f, minimum = %f", fraction, fMinFractionShared));
334  //if(fJetTaggingType==kClosest) Printf("Fraction = %f, minimum = %f", fraction, fMinFractionShared);
336  continue;
337  fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
338  fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
339  if(ptJet2>0.) fh2PtJet2VsRelPt[fCentBin]->Fill(ptJet2,(ptJet1-ptJet2)/ptJet2);
340 
341  Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
342  if(dPhi>TMath::Pi())
343  dPhi -= TMath::TwoPi();
344  if(dPhi<(-1.*TMath::Pi()))
345  dPhi += TMath::TwoPi();
346 
347  fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
348  fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,jet1->DeltaR(jet2));
349  }
350  fNAccJets->Fill(count);
351  return kTRUE;
352  }
353 
355  for(auto j : c.all()){
356  switch(fJetTaggingType){
357  case kClosest:
358  j->ResetMatching();
359  break;
360  case kTag:
361  j->SetTaggedJet(nullptr);
362  j->SetTagStatus(-1);
363  };
364  }
365  }
366 
368  const Int_t kNacceptedBase = contBase.GetNAcceptedJets(),
369  kNacceptedTag = contTag.GetNAcceptedJets();
370  if(!(kNacceptedBase && kNacceptedTag)) return false;
371 
372  // Build kd-trees
373  TArrayD etaBase(kNacceptedBase), phiBase(kNacceptedBase),
374  etaTag(kNacceptedTag), phiTag(kNacceptedTag);
375  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
376 
377  int countBase(0), countTag(0);
378  for(auto jb : contBase.accepted()) {
379  etaBase[countBase] = jb->Eta();
380  phiBase[countBase] = jb->Phi();
381  jetsBase[countBase] = jb;
382  countBase++;
383  }
384  for(auto jt : contTag.accepted()) {
385  etaTag[countTag] = jt->Eta();
386  phiTag[countTag] = jt->Phi();
387  jetsTag[countTag] = jt;
388  countTag++;
389  }
390  TKDTreeID treeBase(etaBase.GetSize(), 2, 1), treeTag(etaTag.GetSize(), 2, 1);
391  treeBase.SetData(0, etaBase.GetArray());
392  treeBase.SetData(1, phiBase.GetArray());
393  treeBase.Build();
394  treeTag.SetData(0, etaTag.GetArray());
395  treeTag.SetData(1, phiTag.GetArray());
396  treeTag.Build();
397 
398  TArrayI faMatchIndexTag(kNacceptedBase), faMatchIndexBase(kNacceptedTag);
399  faMatchIndexBase.Reset(-1);
400  faMatchIndexTag.Reset(-1);
401 
402  // find the closest distance to the full jet
403  countBase = 0;
404  for(auto j : contBase.accepted()) {
405  Double_t point[2] = {j->Eta(), j->Phi()};
406  Int_t index(-1); Double_t distance(-1);
407  treeTag.FindNearestNeighbors(point, 1, &index, &distance);
408  // test whether indices are matching:
409  if(index >= 0 && distance < maxDist){
410  AliDebugStream(1) << "Found closest tag jet for " << countBase << " with match index " << index << " and distance " << distance << std::endl;
411  faMatchIndexTag[countBase]=index;
412  } else {
413  AliDebugStream(1) << "Not found closest tag jet for " << countBase << ", distance to closest " << distance << std::endl;
414  }
415 
416 #ifdef JETTAGGERFAST_TEST
417  if(index>-1){
418  Double_t distanceTest(-1);
419  distanceTest = TMath::Sqrt(TMath::Power(etaTag[index] - j->Eta(), 2) + TMath::Power(phiTag[index] - j->Phi(), 2));
420  if(TMath::Abs(distanceTest - distance) > DBL_EPSILON){
421  AliDebugStream(1) << "Mismatch in distance from tag jet with index from tree: " << distanceTest << ", distance from tree " << distance << std::endl;
422  fIndexErrorRateBase->Fill(1);
423  }
424  }
425 #endif
426 
427  countBase++;
428  }
429 
430  // other way around
431  countTag = 0;
432  for(auto j : contTag.accepted()){
433  Double_t point[2] = {j->Eta(), j->Phi()};
434  Int_t index(-1); Double_t distance(-1);
435  treeBase.FindNearestNeighbors(point, 1, &index, &distance);
436  if(index >= 0 && distance < maxDist){
437  AliDebugStream(1) << "Found closest base jet for " << countBase << " with match index " << index << " and distance " << distance << std::endl;
438  faMatchIndexBase[countTag]=index;
439  } else {
440  AliDebugStream(1) << "Not found closest tag jet for " << countBase << ", distance to closest " << distance << std::endl;
441  }
442 
443 #ifdef JETTAGGERFAST_TEST
444  if(index>-1){
445  Double_t distanceTest(-1);
446  distanceTest = TMath::Sqrt(TMath::Power(etaBase[index] - j->Eta(), 2) + TMath::Power(phiBase[index] - j->Phi(), 2));
447  if(TMath::Abs(distanceTest - distance) > DBL_EPSILON){
448  AliDebugStream(1) << "Mismatch in distance from base jet with index from tree: " << distanceTest << ", distance from tree " << distance << std::endl;
449  fIndexErrorRateTag->Fill(1);
450  }
451  }
452 #endif
453 
454  countTag++;
455  }
456 
457  // check for "true" correlations
458  // these are pairs where the base jet is the closest to the tag jet and vice versa
459  // As the lists are linear a loop over the outer base jet is sufficient.
460  AliDebugStream(1) << "Starting true jet loop: nbase(" << kNacceptedBase << "), ntag(" << kNacceptedTag << ")\n";
461  for(int ibase = 0; ibase < kNacceptedBase; ibase++) {
462  AliDebugStream(2) << "base jet " << ibase << ": match index in tag jet container " << faMatchIndexTag[ibase] << "\n";
463  if(faMatchIndexTag[ibase] > -1){
464  AliDebugStream(2) << "tag jet " << faMatchIndexTag[ibase] << ": matched base jet " << faMatchIndexBase[faMatchIndexTag[ibase]] << "\n";
465  }
466  if(faMatchIndexTag[ibase] > -1 && faMatchIndexBase[faMatchIndexTag[ibase]] == ibase) {
467  AliDebugStream(2) << "found a true match \n";
468  AliEmcalJet *jetBase = jetsBase[ibase],
469  *jetTag = jetsTag[faMatchIndexTag[ibase]];
470  if(jetBase && jetTag) {
471 #ifdef JETTAGGERFAST_TEST
472  if(TMath::Abs(etaBase[ibase] - jetBase->Eta()) > DBL_EPSILON || TMath::Abs(phiBase[ibase] - jetBase->Phi()) > DBL_EPSILON){
473  AliErrorStream() << "Selected incorrect base jet for tagging : eta test(" << jetBase->Eta() << ")/true(" << etaBase[ibase]
474  << "), phi test(" << jetBase->Phi() << ")/true(" << phiBase[ibase] << ")\n";
475  fContainerErrorRateBase->Fill(1);
476  }
477  if(TMath::Abs(etaTag[faMatchIndexTag[ibase]] - jetTag->Eta()) > DBL_EPSILON || TMath::Abs(phiTag[faMatchIndexTag[ibase]] - jetTag->Phi()) > DBL_EPSILON){
478  AliErrorStream() << "Selected incorrect tag jet for tagging : eta test(" << jetTag->Eta() << ")/true(" << etaTag[faMatchIndexTag[ibase]]
479  << "), phi test(" << jetTag->Phi() << ")/true(" << phiTag[faMatchIndexTag[ibase]] << ")\n";
480  fContainerErrorRateTag->Fill(1);
481  }
482 #endif
483  // Test if the position of the jets correp
484  Double_t dR = jetBase->DeltaR(jetTag);
485  switch(fJetTaggingType){
486  case kTag:
487  jetBase->SetTaggedJet(jetTag);
488  jetBase->SetTagStatus(1);
489 
490  jetTag->SetTaggedJet(jetBase);
491  jetTag->SetTagStatus(1);
492  break;
493  case kClosest:
494  jetBase->SetClosestJet(jetTag,dR);
495  jetTag->SetClosestJet(jetBase,dR);
496  break;
497  };
498  }
499  }
500  }
501  return kTRUE;
502  }
503 
505  return GetDeltaPhi(jet1->Phi(),jet2->Phi());
506  }
507 
509  Double_t dPhi = phi1-phi2;
510  if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
511  if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
512 
513  return dPhi;
514  }
515 
517  const char * njetsTag,
518  const Double_t R,
519  const char * nrhoBase,
520  const char * nrhoTag,
521  const char * ntracks,
522  const char * nclusters,
523  const char * type,
524  const char * CentEst,
525  Int_t pSel,
526  const char * trigClass){
527 
528  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
529  if (!mgr) {
530  std::cerr << "E-AddEmcalJetTaggerTaskFast: No analysis manager found.\n";
531  return nullptr;
532  }
533  // Check the analysis type using the event handlers connected to the analysis manager.
534  //==============================================================================
535  if (!mgr->GetInputEventHandler())
536  {
537  std::cerr << "E-AddEmcalJetTaggerTaskFast: This task requires an input event handler\n";
538  return NULL;
539  }
540 
541  TString wagonName = Form("JetTagger_%s_%s_TC%s",njetsBase,njetsTag,trigClass);
542 
543  //Configure jet tagger task
545 
546  task->SetNCentBins(4);
547  //task->SetVzRange(-10.,10.);
548 
549  AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
550  AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
551 
552  task->SetJetContainerBase(0);
553  task->SetJetContainerTag(1);
554 
555  TString strType(type);
556  AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
557  if(jetContBase) {
558  jetContBase->SetRhoName(nrhoBase);
559  jetContBase->ConnectParticleContainer(trackCont);
560  jetContBase->ConnectClusterContainer(clusterCont);
561  jetContBase->SetMaxTrackPt(10000.);
562  }
563 
564  AliJetContainer *jetContTag = task->AddJetContainer(njetsTag,"TPC",R);
565  if(jetContTag) {
566  jetContTag->SetRhoName(nrhoTag);
567  jetContTag->ConnectParticleContainer(trackCont);
568  jetContTag->ConnectClusterContainer(clusterCont);
569  jetContTag->SetMaxTrackPt(10000.);
570  }
571  for(Int_t i=0; i<2; i++) {
572  task->SetPercAreaCut(0.6, i); //keep?
573  }
574  task->SetCentralityEstimator(CentEst);
575  task->SelectCollisionCandidates(pSel);
576  task->SetUseAliAnaUtils(kFALSE);
577 
578  mgr->AddTask(task);
579 
580  //Connnect input
581  mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
582 
583  //Connect output
584  TString contName(wagonName);
585  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
586  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
587  mgr->ConnectOutput(task,1,coutput1);
588 
589  return task;
590  }
591 
592 
593  }
594 }
595 
void SetTagStatus(Int_t i)
Definition: AliEmcalJet.h:350
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:351
Definition: External.C:260
Definition: External.C:236
Double_t GetJetEtaMin() const
void SetTaggedJet(AliEmcalJet *j)
Definition: AliEmcalJet.h:349
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliJetContainer * GetJetContainer(Int_t i=0) const
Int_t GetTagStatus() const
Definition: AliEmcalJet.h:352
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:337
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)
Enable general histograms.
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.
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.
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