AliPhysics  9fe175b (9fe175b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetMatching.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 //
17 // General task to match two sets of jets
18 //
19 // This task takes two TClonesArray's as input:
20 // [1] fSourceJets - e.g. pythia jets
21 // [2] fTargetJets - e.g. a samle containing pythia jets embedded in a pbpb event
22 // the task will try to match jets from the source array to the target array, and
23 // save the found TARGET jets in a new array 'fMatchedJets'
24 // the purpose of this task is twofold
25 // [1] matching for embedding studies
26 // [2] matching to create a detector response function
27 //
28 // matching is done in three steps
29 // [1] geometric matching, where jets are matched by R = sqrt(dphi*dphi-deta*deta) or directly via dphi and deta
30 // [2] optional injection / bijection check
31 // in this check, it is made sure that fSourceJets -> fMatchedJets is either an injective non-surjective
32 // or bijective map, depending on the matching resolution which is chosen for step [1]
33 // so that each source jet has a unique target and vice-versa.
34 // if R (or dphi, deta) are proportional to, or larger than, the jet radius, matching tends to be
35 // bijective (each source has a target), if R is chosen to be very small, source jets might be lost
36 // as no target can be found.
37 // the mapping is done in such a way that each target is matched to its closest source and each source
38 // is mapped to its closest target
39 // [3] constituent matching
40 // - how many constituents of the source jet are present in the matched jet?
41 // a cut on the constituent fraction can be performed (not recommended)
42 // - how much of the original pt is recovered in the matched jet?
43 // a cut on the fraction of recovered / original pt can be performed (recommended)
44 //
45 // detector response
46 // to obtain a detector respose function, supply
47 // [1] fSourceJets - particle level jets
48 // [2] fTargetJets - detector level jets
49 //
50 // Author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
51 // rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
52 
53 // root includes
54 #include <TClonesArray.h>
55 #include <TChain.h>
56 #include <TMath.h>
57 #include <TF1.h>
58 #include <TF2.h>
59 #include <TH3.h>
60 #include <TH1F.h>
61 #include <TH2F.h>
62 #include <TH3F.h>
63 #include <TProfile.h>
64 #include <TFile.h>
65 #include <TTree.h>
66 #include <TKey.h>
67 #include <TSystem.h>
68 // aliroot includes
69 #include <AliAnalysisTask.h>
70 #include <AliAnalysisManager.h>
71 #include <AliLog.h>
72 #include <AliVEvent.h>
73 #include <AliVParticle.h>
74 // emcal jet framework includes
75 #include <AliEmcalJet.h>
77 #include <AliLocalRhoParameter.h>
78 
80 using namespace std;
81 
83 
86  // default constructor
88 }
89 //_____________________________________________________________________________
92  // constructor
94  DefineInput(0, TChain::Class());
95  DefineOutput(1, TList::Class());
96 }
97 //_____________________________________________________________________________
99 {
100  // destructor
101  if(fOutputList) delete fOutputList;
102 }
103 //_____________________________________________________________________________
105 {
106  // initialize the anaysis
107  #ifdef DEBUGTASK
108  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
109  #endif
110  // get the stand alone jets from the input event
111  fSourceJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSourceJetsName.Data()));
112  if(!fSourceJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
113  fTargetJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTargetJetsName.Data()));
114  if(!fTargetJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
115  // append the list of matched jets to the event
116  fMatchedJets->Delete();
117  if (!(InputEvent()->FindListObject(fMatchedJetsName))) InputEvent()->AddObject(fMatchedJets);
118  else AliFatal(Form("%s: Object with name %s already in event! Aborting", GetName(), fMatchedJetsName.Data()));
120  switch (fSourceBKG) {
121  case kSourceLocalRho : {
122  fSourceRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fSourceRhoName));
123  if(!fSourceRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fSourceRhoName.Data()));
124  } break;
125  default : break;
126  }
127  switch (fTargetBKG) {
128  case kTargetLocalRho : {
129  fTargetRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fTargetRhoName));
130  if(!fTargetRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fTargetRhoName.Data()));
131  } break;
132  default : break;
133  }
134  AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
136 }
137 //_____________________________________________________________________________
139 {
140  // create output objects
141  #ifdef DEBUGTASK
142  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
143  #endif
144  // add the matched jets array to the event
145  fMatchedJets = new TClonesArray("AliEmcalJet");
146  fMatchedJets->SetName(fMatchedJetsName);
147  fOutputList = new TList();
148  fOutputList->SetOwner(kTRUE);
149  // add analysis histograms
150  fHistUnsortedCorrelation = BookTH1F("fHistUnsortedCorrelation", "#Delta #varphi unsorted", 50, 0, TMath::Pi());
151  fHistMatchedCorrelation = BookTH1F("fHistMatchedCorrelation", "#Delta #varphi matched", 50, 0, TMath::Pi());
152  fHistSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
153  fHistMatchedSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistMatchedParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistMatchedSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
154  fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
155  fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
156  fHistSourceMatchedJetPt = (fDoDetectorResponse) ? BookTH2F("fHistDetectorResponse", "particle level jet p_{t}^{gen} [GeV/c]", "detector level jet p_{t}^{rec} [GeV/c]", 300, -50, 250, 300, -50, 250) : BookTH2F("fHistSourceMatchedJetPt", "source jet p_{t} [GeV/c]", "matched jet p_{t} [GeV/c]", 300, -50, 250, 300, -50, 250);
157  if(fDoDetectorResponse) {
158  fHistDetectorResponseProb = BookTH2F("fHistDetectorResponseProb", "(p_{t}^{det} - p_{t}^{part})/p_{t}^{part}", "p_{t}^{part}", 100, -1.5, 1., 20, 0, 200);
159  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
160  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
161  fOutputList->Add(fh1Xsec);
162  fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
163  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
164  fOutputList->Add(fh1Trials);
165  fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
166  fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
168  }
169  fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents source jets", 50, 0, 50);
170  fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents target jets", 50, 0, 50);
171  fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents matched jets", 50, 0, 50);
172  if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
173  fProfFracPtMatched = new TProfile("fProfFracPtMatched", "recovered target p_{T} / source p_{T}", 15, -50, 250);
175  fProfFracNoMatched = new TProfile("fProfFracNoMatched", "recovered target constituents / source constituents", 15, -50, 250);
177  }
178  // the analysis summary histo which stores all the analysis flags is always written to file
179  fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
180  fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
181  fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
182  fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
183  fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
185  fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
186  fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
187  fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
188  fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
189  fOutputList->Add(fProfQA);
190  fProfFracPtJets = new TProfile("fProfFracPtJets", "source p_{T} / target p_{T}", 15, -50, 250);
192  fProfFracNoJets = new TProfile("fProfFracNoJets", "source constituents / target constituents", 15, -50, 250);
194  switch (fMatchingScheme) {
195  case kDiJet : {
196  fHistDiJet = BookTH3F("fHistDiJet", "matched di-jet #varphi", "matched di-jet #eta", "leading jet p_{T} (GeV/c)", 100, 0., TMath::TwoPi(), 100, -5., 5., 100, 0, 200);
197  fHistDiJetLeadingJet = BookTH3F("fHistDiJetLeadingJet", "leading jet #varphi", "leadingd jet #eta", "leading jet p_{T} (GeV/c)", 100, 0., TMath::TwoPi(), 100, -5., 5., 100, 0, 200);
198  fHistDiJetDPhi = BookTH2F("fHistDiJetDPhi", "leading jet #varphi - (matched jet #varphi - #pi)", "leading jet p_{T} (GeV/c)", 100, -1.*TMath::Pi(), TMath::Pi(), 100, 0, 200);
199  fHistDiJetDPt = BookTH2F("fHistDiJetDPt", "leading jet p_{T} - sub leading jet p_{T} (GeV/c)", "leading jet p_{T} (GeV/c)", 100, -25, 25, 100, 0, 200);
200  } break;
201  default : break;
202  }
203  fOutputList->Sort();
204  PostData(1, fOutputList);
205 }
206 //_____________________________________________________________________________
207 TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
208 {
209  // book a TH1F and connect it to the output container
210  #ifdef DEBUGTASK
211  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
212  #endif
213  if(append && !fOutputList) return 0x0;
214  TString title(name);
215  title += Form(";%s;[counts]", x);
216  TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
217  histogram->Sumw2();
218  if(append) fOutputList->Add(histogram);
219  return histogram;
220 }
221 //_____________________________________________________________________________
222 TH2F* AliAnalysisTaskJetMatching::BookTH2F(const char* name, const char* x, const char*y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Bool_t append)
223 {
224  // book a TH2F and connect it to the output container
225  #ifdef DEBUGTASK
226  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
227  #endif
228  if(append && !fOutputList) return 0x0;
229  TString title(name);
230  title += Form(";%s;%s", x, y);
231  TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
232  histogram->Sumw2();
233  if(append) fOutputList->Add(histogram);
234  return histogram;
235 }
236 //_____________________________________________________________________________
237 TH3F* AliAnalysisTaskJetMatching::BookTH3F(const char* name, const char* x, const char* y, const char* z, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t binsz, Double_t minz, Double_t maxz, Bool_t append)
238 {
239  // book a TH2F and connect it to the output container
240  #ifdef DEBUGTASK
241  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
242  #endif
243  if(append && !fOutputList) return 0x0;
244  TString title(name);
245  title += Form(";%s;%s;%s", x, y, z);
246  TH3F* histogram = new TH3F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy, binsz, minz, maxz);
247  histogram->Sumw2();
248  if(append) fOutputList->Add(histogram);
249  return histogram;
250 }
251 //_____________________________________________________________________________
253 {
254  // for each file get the number of trials and pythia cross section
255  // see $ALICE_PHYSICS/PWGJE/AliAnalysisTaskJetProperties.cxx
256  // $ALICE_PHYSICS/PWG/Tools/AliAnalysisHelperJetTasks.cxx
257  // this function is implenented here temporarily to avoid introducing a dependency
258  // later on this could just be a call to a static helper function
259  #ifdef DEBUGTASK
260  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
261  #endif
262  Float_t xsection(0), ftrials(1);
263  fAvgTrials = ftrials;
264  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
265  if(tree) {
266  TFile *curfile = tree->GetCurrentFile();
267  if (!curfile) return kFALSE;
268  TString file(curfile->GetName());
269  if(file.Contains("root_archive.zip#")) file.Replace(file.Index("#",1,file.Index("root_archive",12,0,TString::kExact),TString::kExact)+1,file.Index(".root",5,TString::kExact)-file.Index("root_archive",12,0,TString::kExact),"");
270  else file.ReplaceAll(gSystem->BaseName(file.Data()),"");
271  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
272  if(!fxsec) {
273  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
274  if(fxsec) {
275  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
276  if(key) {
277  TList *list = dynamic_cast<TList*>(key->ReadObj());
278  if(list) {
279  xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
280  ftrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
281  }
282  }
283  fxsec->Close();
284  }
285  } else {
286  TTree *xtree = (TTree*)fxsec->Get("Xsection");
287  if(xtree) {
288  UInt_t _ftrials = 0;
289  Double_t _xsection = 0;
290  xtree->SetBranchAddress("xsection",&_xsection);
291  xtree->SetBranchAddress("ntrials",&_ftrials);
292  xtree->GetEntry(0);
293  ftrials = _ftrials;
294  xsection = _xsection;
295  }
296  fxsec->Close();
297  }
298  fh1Xsec->Fill("<#sigma>",xsection);
299  Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
300  if(ftrials >= nEntries && nEntries > 0.) fAvgTrials = ftrials/nEntries;
301  fh1Trials->Fill("#sum{ntrials}",ftrials);
302  }
303  return kTRUE;
304 }
305 //_____________________________________________________________________________
307 {
308  // execute once for each event
309  #ifdef DEBUGTASK
310  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
311  #endif
312  if(!(InputEvent() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
313  if(fh1AvgTrials) fh1AvgTrials->Fill("#sum{avg ntrials}", fAvgTrials);
314  // step one: do geometric matching
315  switch (fMatchingScheme) {
316  // FIXME having separate dedicated functions is not necessary and historical
317  // cluttered code - should be cleaned up and merged into one
318  case kGeoEtaPhi : {
320  // break if no jet was matched
321  if(!fMatchedJetContainer[1][0]) return kTRUE;
322  } break;
323  case kGeoR : {
325  // break if no jet was matched
326  if(!fMatchedJetContainer[1][0]) return kTRUE;
327  } break;
328  case kDiJet : {
329  DoDiJetMatching();
330  } break;
331  default : break;
332  }
333  // optional step two: get a bijection (avoid duplicate matches)
335  // optional step three: match constituents within matched jets
337  // stream data to output
338  PostMatchedJets();
340  #ifdef DEBUGTASK
341  PrintInfo();
342  #endif
343  PostData(1, fOutputList);
344  return kTRUE;
345 }
346 //_____________________________________________________________________________
348 {
349  // do geometric matching based on eta phi distance between jet centers
350  #ifdef DEBUGTASK
351  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
352  #endif
353  fNoMatchedJets = 0; // reset the matched jet counter
354  Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
355  for(Int_t i(0); i < iSource; i++) {
356  AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
357  if(!PassesCuts(sourceJet, 0)) continue;
358  for(Int_t j(0); j < iTarget; j++) {
359  AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
360  if(!PassesCuts(targetJet, 1)) continue;
361  if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
362  if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
363  Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
364  if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
365  if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
366  if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) { // accept the jets as matching
367  Bool_t isBestMatch(kTRUE);
368  if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
369  #ifdef DEBUGTASK
370  printf(" > Entering first bijection test \n");
371  #endif
372  for(Int_t k(i); k < iSource; k++) {
373  AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
374  if(PassesCuts(candidateSourceJet, 0)) continue;
375  #ifdef DEBUGTASK
376  printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
377  #endif
378  if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
379  isBestMatch = kFALSE;
380  break;
381  }
382  }
383  #ifdef DEBUGTASK
384  (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
385  #endif
386  }
387  if(isBestMatch) {
388  fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
389  fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
390  fNoMatchedJets++;
391  }
392  if(fNoMatchedJets > 199) { // maximum to keep the cache at reasonable size
393  AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
394  return;
395  }
396  }
397  }
398  }
399  }
400 }
401 //_____________________________________________________________________________
403 {
404  // do geometric matching based on shortest path between jet centers
405  #ifdef DEBUGTASK
406  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
407  #endif
408  fNoMatchedJets = 0; // reset the matched jet counter
409  Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
410  for(Int_t i(0); i < iSource; i++) {
411  AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
412  if(!PassesCuts(sourceJet, 0)) continue;
413  for(Int_t j(0); j < iTarget; j++) {
414  AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
415  if(!PassesCuts(targetJet, 1)) continue;
416  else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
417  if(GetR(sourceJet, targetJet) <= fMatchR) {
418  Bool_t isBestMatch(kTRUE);
419  if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
420  #ifdef DEBUGTASK
421  printf(" > Entering first bijection test \n");
422  #endif
423  for(Int_t k(i); k < iSource; k++) {
424  AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
425  if(!PassesCuts(candidateSourceJet, 0)) continue;
426  #ifdef DEBUGTASK
427  printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
428  #endif
429  if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
430  isBestMatch = kFALSE;
431  break;
432  }
433  }
434  #ifdef DEBUGTASK
435  (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
436  #endif
437  }
438  if(isBestMatch) {
439  fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
440  fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
441  fNoMatchedJets++;
442  }
443  if(fNoMatchedJets > 99) {
444  AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
445  return;
446  }
447  }
448  }
449  }
450 }
451 //_____________________________________________________________________________
453 {
454  // match dijets. this is in a sense a 'special' mode of the task as both jet supplied jet arrays
455  // (target and source) are the same jet collection, matching will be performed on distribution in
456  // azimuth
457  // no ouptut array is produced, dedicated dijet histo's are filled in this function instead
458  #ifdef DEBUGTASK
459  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
460  #endif
461  // get the leading jet in a given acceptance (TPC is default)
462  Int_t leadingJetIndex(-1), subLeadingJetIndex(-1);
463  // retrieve the leading jet, leadingJetIndex points to the leading jet
464  AliEmcalJet* leadingJet(GetLeadingJet(fSourceJets, leadingJetIndex));
465  if(!leadingJet) return;
466  // fill phi and eta of leading jet (should always be in selected acceptance)
467  fHistDiJetLeadingJet->Fill(leadingJet->Phi(), leadingJet->Eta(), leadingJet->Pt());
468  Double_t sourcePhi(leadingJet->Phi()), targetPhi(-1);
469  // get the sub-leading jet - faster when leading jet is also provided
470  AliEmcalJet* subLeadingJet(GetSubLeadingJet(fSourceJets, leadingJetIndex, subLeadingJetIndex));
471  if(!subLeadingJet) return;
472  else { // check if the sub-leading jet is actually the away-side jet
473  targetPhi = subLeadingJet->Phi() + TMath::Pi();
474  // rotate jets to common phase
475  if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
476  if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
477  if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi+TMath::TwoPi())) targetPhi+=TMath::TwoPi();
478  if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi-TMath::TwoPi())) targetPhi-=TMath::TwoPi();
479  if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {
480  fHistDiJet->Fill(subLeadingJet->Phi(), subLeadingJet->Eta(), leadingJet->Pt());
481  fHistDiJetDPhi->Fill(sourcePhi-targetPhi, leadingJet->Pt());
482  fHistDiJetDPt->Fill(leadingJet->Pt() - subLeadingJet->Pt(), leadingJet->Pt());
483  }
484  }
485 }
486 //_____________________________________________________________________________
488 {
489  // match constituents within matched jets
490  #ifdef DEBUGTASK
491  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
492  #endif
493  if(!fTracks) {
494  AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
495  return; // coverity ...
496  }
497  for(Int_t i(0); i < fNoMatchedJets; i++) {
498  AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
499  AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
500  if(sourceJet && targetJet) { // duplicate check: slot migth be NULL
501  Double_t targetPt(0);
502  Int_t iSJ(sourceJet->GetNumberOfTracks());
503  Int_t iTJ(targetJet->GetNumberOfTracks());
504  Int_t overlap(0), alreadyFound(0);
505  for(Int_t j(0); j < iSJ; j++) {
506  alreadyFound = 0;
507  Int_t idSource((Int_t)sourceJet->TrackAt(j));
508  for(Int_t k(0); k < iTJ; k++) { // compare all target tracks to the source track
509  if(idSource == targetJet->TrackAt(k) && alreadyFound == 0) {
510  overlap++;
511  alreadyFound++; // avoid possible duplicate matching
512  AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
513  if(vp) targetPt += vp->Pt();
514  continue;
515  }
516  }
517  }
518  if((float)overlap/(float)iSJ < fMinFracRecoveredConstituents || targetPt/sourceJet->Pt() < fMinFracRecoveredConstituentPt) {
519  #ifdef DEBUGTASK
520  printf(" \n > Purging jet, recovered constituents ratio %i / %i = %.2f < or pt ratio %.2f / %.2f = %.2f < %.2f", overlap, iSJ, (float)overlap/(float)iSJ, targetPt, sourceJet->Pt(), targetPt/sourceJet->Pt(), fMinFracRecoveredConstituentPt);
521  #endif
522  fMatchedJetContainer[i][0] = 0x0;
523  fMatchedJetContainer[i][1] = 0x0;
524  continue;
525  }
526  if(sourceJet->Pt() > 0) {
527  Double_t sourceRho(0), targetRho(0);
528  if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(sourceJet->Phi(), fSourceRadius)*sourceJet->Area();
529  if(fTargetRho) targetRho = fTargetRho->GetLocalVal(targetJet->Phi(), fTargetRadius)*targetJet->Area();
530  fProfFracPtMatched->Fill(sourceJet->Pt()-sourceRho, (targetPt-targetRho) / (sourceJet->Pt()-sourceRho));
531  fProfFracPtJets->Fill(sourceJet->Pt()-sourceRho, (targetJet->Pt()-targetRho) / (sourceJet->Pt()-sourceRho));
532  fProfFracNoMatched->Fill(sourceJet->Pt()-sourceRho, (double)overlap / (double)sourceJet->GetNumberOfTracks());
533  fProfFracNoJets->Fill(sourceJet->Pt()-sourceRho, (double)targetJet->GetNumberOfTracks() / (double)sourceJet->GetNumberOfTracks());
534  }
535  #ifdef DEBUGTASK
536  if(fDebug > 0) {
537  printf("\n > Jet A: %i const\t", iSJ);
538  printf(" > Jet B %i const\t", iTJ);
539  printf(" -> OVERLAP: %i tracks <- \n", overlap);
540  }
541  #endif
542  }
543  }
544 }
545 //_____________________________________________________________________________
547 {
548  // bijection of source and matched jets, based on closest distance between jets
549  #ifdef DEBUGTASK
550  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
551  #endif
552  for(Int_t i(0); i < fNoMatchedJets; i++) {
553  for(Int_t j(i+1); j < fNoMatchedJets; j++) {
554  if(fMatchedJetContainer[i][0] == fMatchedJetContainer[j][0]) {
555  // found source with two targets, now see which target is closer to the source
556  if(!(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1] && fMatchedJetContainer[j][0] && fMatchedJetContainer[j][1] )) continue;
557  Double_t rA(GetR(fMatchedJetContainer[i][0], fMatchedJetContainer[i][1])); // distance between connection A = i
558  Double_t rB(GetR(fMatchedJetContainer[j][0], fMatchedJetContainer[j][1])); // distance between connection B = j
559  if (rB > rA) { // jet two is far away, clear it from both target and source list
560  fMatchedJetContainer[j][0] = 0x0;
561  fMatchedJetContainer[j][1] = 0x0;
562  } else { // jet one is far away, clear it from both target and source list
563  fMatchedJetContainer[i][0] = fMatchedJetContainer[j][0]; // switch to avoid breaking loop
565  fMatchedJetContainer[j][0] = 0x0; // clear duplicate jet from cache
566  fMatchedJetContainer[j][1] = 0x0;
567  }
568  #ifdef DEBUGTASK
569  printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
570  #endif
571  }
572  }
573  }
574 }
575 //_____________________________________________________________________________
577 {
578  // post matched jets
579  #ifdef DEBUGTASK
580  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
581  #endif
582  fMatchedJets->Delete(); // should be a NULL operation, but added just in case
583  for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
584  if(fMatchedJetContainer[i][1]) { // duplicate jets will have NULL value here and are skipped
585  new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
586  p++;
587  }
588  }
589 }
590 //_____________________________________________________________________________
592 {
593  // fill the analysis summary histrogram, saves all relevant analysis settigns
594  #ifdef DEBUGTASK
595  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
596  #endif
597  fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
598  fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
599  fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
600  fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
601  fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
602  fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
603  fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
604  fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
605  fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
606  fHistAnalysisSummary->SetBinContent(5, fMatchEta);
607  fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
608  fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
609  fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
610  fHistAnalysisSummary->SetBinContent(7, fMatchR);
611  fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
612  fHistAnalysisSummary->SetBinContent(8, 1.);
613 }
614 //_____________________________________________________________________________
616 {
617  // fill matched jet histograms
618  #ifdef DEBUGTASK
619  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
620  #endif
621  for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
622  AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
623  if(!source) continue;
624  else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
625  Double_t sourceRho(0), targetRho(0);
626  if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
627  fHistSourceJetPt->Fill(source->Pt()-sourceRho);
629  for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
630  AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
631  if(target) {
632  if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
633  if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
634  fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));
635  fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
636  fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
637  fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
638  if(j==0) {
639  fHistTargetJetPt->Fill(target->Pt()-targetRho);
641  }
642  }
643  }
644  }
645  for(Int_t i(0); i < fNoMatchedJets; i++) {
646  if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
647  Double_t sourceRho(0), targetRho(0);
651  fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
652  fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
653  fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
654  fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
655  fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
656  fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
657 
658  fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
659  if(fDoDetectorResponse) {
660  fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
661  fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
662  fHistDetectorResponseProb->Fill((fMatchedJetContainer[i][1]->Pt()-fMatchedJetContainer[i][0]->Pt())/fMatchedJetContainer[i][0]->Pt(), fMatchedJetContainer[i][0]->Pt());
663  }
664  }
665  }
666 }
667 //_____________________________________________________________________________
668 AliEmcalJet* AliAnalysisTaskJetMatching::GetLeadingJet(TClonesArray* source, Int_t &leadingJetIndex, Double_t etaMin, Double_t etaMax)
669 {
670  // return the leading jet within giiven acceptance
671  Int_t iJets(source->GetEntriesFast());
672  Double_t pt(0);
673  AliEmcalJet* leadingJet(0x0);
674  for(Int_t i(0); i < iJets; i++) {
675  AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
676  if(jet->Eta() < etaMin || jet->Eta() > etaMax) continue;
677  if(jet->Pt() > pt) {
678  leadingJet = jet;
679  pt = leadingJet->Pt();
680  leadingJetIndex = i;
681  }
682  }
683  return leadingJet;
684 }
685 //_____________________________________________________________________________
686 AliEmcalJet* AliAnalysisTaskJetMatching::GetSubLeadingJet(TClonesArray* source, Int_t leadingJetIndex, Int_t &subLeadingJetIndex)
687 {
688  // return the sub-leading jet within given acceptance
689  // same as GetLeadingJet() but skips the leading jet (so returned jet is
690  // sub-leading by design)
691  // there is no eta requirement on the location of the sub-leading jet
692  Int_t iJets(source->GetEntriesFast());
693  // if the leading jet isn't given, retrieve it
694  if(leadingJetIndex < 0) GetLeadingJet(source, leadingJetIndex);
695  AliEmcalJet *leadingJet(0x0), *prevLeadingJet(static_cast<AliEmcalJet*>(source->At(leadingJetIndex)));
696  if(!prevLeadingJet) return 0x0;
697  Double_t pt(0), leadingPt(prevLeadingJet->Pt());
698  for(Int_t i(0); i < iJets; i++) {
699  if(i == leadingJetIndex) continue;
700  AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
701  // check if jet is actually sub-leading
702  if(jet->Pt() > pt && jet->Pt() <= leadingPt) {
703  leadingJet = jet;
704  pt = leadingJet->Pt();
705  subLeadingJetIndex = i;
706  }
707  }
708  return leadingJet;
709 }
710 //_____________________________________________________________________________
712 {
713  // print some info
714  printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
715  printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
716  printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
717 }
718 //_____________________________________________________________________________
720 {
721  // terminate
722 }
723 //_____________________________________________________________________________
TH1F * fHistTargetJetPt
pt of matched source jets
TH3F * fHistDiJet
no of consstituents fraction jet / jet
TH3F * fHistDiJetLeadingJet
matched dijet eta, phi
Bool_t PassesCuts(const AliVTrack *track) const
Double_t Area() const
Definition: AliEmcalJet.h:114
TH1F * fHistMatchedJetPt
pt of target jets
const char * title
Definition: MakeQAPdf.C:26
ClassImp(AliAnalysisTaskJetMatching) AliAnalysisTaskJetMatching
Double_t Eta() const
Definition: AliEmcalJet.h:105
TH2F * fHistSourceMatchedJetPt
pt of matched jets
TString fTargetJetsName
array with target jets
Double_t Phi() const
Definition: AliEmcalJet.h:101
TProfile * fProfFracNoMatched
sum pt fraction for matched jets
TSystem * gSystem
Double_t GetLocalVal(Double_t phi, Double_t r, Double_t n) const
TList * list
TH1F * fHistNoConstMatchJet
number of constituents of target jet
TString fMatchedJetsName
final list of matched jets which is added to event
Double_t GetR(AliVParticle *a, AliVParticle *b) const
TH1F * fh1AvgTrials
sum of trails (in Notify)
TH3F * BookTH3F(const char *name, const char *x, const char *y, const char *z, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t binsz, Double_t minz, Double_t maxz, Bool_t append=kTRUE)
AliEmcalJet * fMatchedJetContainer[200][2]
number of matched jets
TH2F * fHistDetectorResponseProb
pt or source vs matched jets
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:124
AliEmcalJet * GetLeadingJet(TClonesArray *source, Int_t &leadingJetIndex, Double_t etaMin=-.7, Double_t etaMax=.7)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:123
TH1F * BookTH1F(const char *name, const char *x, Int_t bins, Double_t min, Double_t max, Bool_t append=kTRUE)
Float_t fMatchEta
pointers to matched jets
TH1F * fHistMatchedSourceJetPt
pt of source jets
Double_t PhaseShift(Double_t x) const
TProfile * fProfFracPtMatched
number of constituents of matched jets
TString fSourceJetsName
array with source jets
TProfile * fProfFracPtJets
sum pt fraction for matched tracks / jet
TProfile * fh1Xsec
average sum of trials (in Run)
TProfile * fProfFracNoJets
no of constituents fraction found / jet
Float_t fAvgTrials
pythia cross section and trials
Double_t Pt() const
Definition: AliEmcalJet.h:93
TProfile * fProfQA
QA spreads of matched jets.
virtual Bool_t IsEventSelected()
virtual void Terminate(Option_t *option)
TClonesArray * fTracks
!tracks
TH2F * fHistDiJetDPhi
leading jet (for dijet) eta, phi
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:144
TH2F * BookTH2F(const char *name, const char *x, const char *y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Bool_t append=kTRUE)
TFile * file
AliEmcalJet * GetSubLeadingJet(TClonesArray *source, Int_t leadingJetIndex, Int_t &subLeadingJetIndex)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Int_t fNoMatchedJets
QA spreads of source and target jets.
TH1F * fHistMatchedCorrelation
dphi correlation of unsorted jets
virtual Bool_t AcceptJet(AliEmcalJet *jet, Int_t c=0)
TH1F * fHistNoConstTargetJet
number of constituents of source jet
TH1F * fHistSourceJetPt
dphi correlation of matched jets