]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Removing printout (Davide)
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetMatching.cxx
... / ...
CommitLineData
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 <TH1F.h>
60#include <TH2F.h>
61#include <TProfile.h>
62#include <TFile.h>
63#include <TTree.h>
64#include <TKey.h>
65#include <TSystem.h>
66// aliroot includes
67#include <AliAnalysisTask.h>
68#include <AliAnalysisManager.h>
69#include <AliLog.h>
70#include <AliVEvent.h>
71#include <AliVParticle.h>
72// emcal jet framework includes
73#include <AliEmcalJet.h>
74#include <AliAnalysisTaskJetMatching.h>
75#include <AliLocalRhoParameter.h>
76
77class AliAnalysisTaskJetMatching;
78using namespace std;
79
80ClassImp(AliAnalysisTaskJetMatching)
81
82AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", kTRUE),
83 fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fSourceRho(0), fSourceRhoName(0), fTargetRho(0), fTargetRhoName(0), fUseScaledRho(0), fSourceRadius(0.3), fTargetRadius(0.3), fMatchingScheme(kGeoEtaPhi), fUseEmcalBaseJetCuts(kFALSE), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistMatchedSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistSourceMatchedJetPt(0), fHistDetectorResponseProb(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fProfFracPtMatched(0), fProfFracPtJets(0), fProfFracNoMatched(0), fProfFracNoJets(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(200), fMatchEta(.3), fMatchPhi(.3), fMatchR(.08), fDoDetectorResponse(kFALSE), fMatchConstituents(kTRUE), fMinFracRecoveredConstituents(.5), fMinFracRecoveredConstituentPt(0.5), fGetBijection(kTRUE), fh1Trials(0x0), fh1AvgTrials(0x0), fh1Xsec(0x0), fAvgTrials(0) {
84 // default constructor
85 ClearMatchedJetsCache();
86}
87//_____________________________________________________________________________
88AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
89 fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fSourceRho(0), fSourceRhoName(0), fTargetRho(0), fTargetRhoName(0), fUseScaledRho(0), fSourceRadius(0.3), fTargetRadius(0.3), fMatchingScheme(kGeoEtaPhi), fUseEmcalBaseJetCuts(kFALSE), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistMatchedSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistSourceMatchedJetPt(0), fHistDetectorResponseProb(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fProfFracPtMatched(0), fProfFracPtJets(0), fProfFracNoMatched(0), fProfFracNoJets(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(200), fMatchEta(.3), fMatchPhi(.3), fMatchR(.08), fDoDetectorResponse(kFALSE), fMatchConstituents(kTRUE), fMinFracRecoveredConstituents(0.5), fMinFracRecoveredConstituentPt(0.5), fGetBijection(kTRUE), fh1Trials(0x0), fh1AvgTrials(0x0), fh1Xsec(0x0), fAvgTrials(0) {
90 // constructor
91 ClearMatchedJetsCache();
92 DefineInput(0, TChain::Class());
93 DefineOutput(1, TList::Class());
94}
95//_____________________________________________________________________________
96AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
97{
98 // destructor
99 if(fOutputList) delete fOutputList;
100}
101//_____________________________________________________________________________
102void AliAnalysisTaskJetMatching::ExecOnce()
103{
104 // initialize the anaysis
105 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
106 // get the stand alone jets from the input event
107 fSourceJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSourceJetsName.Data()));
108 if(!fSourceJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
109 fTargetJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTargetJetsName.Data()));
110 if(!fTargetJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
111 // append the list of matched jets to the event
112 fMatchedJets->Delete();
113 if (!(InputEvent()->FindListObject(fMatchedJetsName))) InputEvent()->AddObject(fMatchedJets);
114 else AliFatal(Form("%s: Object with name %s already in event! Aborting", GetName(), fMatchedJetsName.Data()));
115 FillAnalysisSummaryHistogram();
116 switch (fSourceBKG) {
117 case kSourceLocalRho : {
118 fSourceRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fSourceRhoName));
119 if(!fSourceRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fSourceRhoName.Data()));
120 } break;
121 default : break;
122 }
123 switch (fTargetBKG) {
124 case kTargetLocalRho : {
125 fTargetRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fTargetRhoName));
126 if(!fTargetRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fTargetRhoName.Data()));
127 } break;
128 default : break;
129 }
130 AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
131 if(fDoDetectorResponse) fMatchConstituents = kFALSE;
132}
133//_____________________________________________________________________________
134void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
135{
136 // create output objects
137 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
138 // add the matched jets array to the event
139 fMatchedJets = new TClonesArray("AliEmcalJet");
140 fMatchedJets->SetName(fMatchedJetsName);
141 fOutputList = new TList();
142 fOutputList->SetOwner(kTRUE);
143 // add analysis histograms
144 fHistUnsortedCorrelation = BookTH1F("fHistUnsortedCorrelation", "#Delta #varphi unsorted", 50, 0, TMath::Pi());
145 fHistMatchedCorrelation = BookTH1F("fHistMatchedCorrelation", "#Delta #varphi matched", 50, 0, TMath::Pi());
146 fHistSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
147 fHistMatchedSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistMatchedParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistMatchedSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
148 fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
149 fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
150 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);
151 if(fDoDetectorResponse) {
152 fHistDetectorResponseProb = BookTH2F("fHistDetectorResponseProb", "(p_{t}^{det} - p_{t}^{part})/p_{t}^{part}", "p_{t}^{part}", 100, -1.5, 1., 20, 0, 200);
153 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
154 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
155 fOutputList->Add(fh1Xsec);
156 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
157 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
158 fOutputList->Add(fh1Trials);
159 fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
160 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
161 fOutputList->Add(fh1AvgTrials);
162 }
163 fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents source jets", 50, 0, 50);
164 fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents target jets", 50, 0, 50);
165 fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents matched jets", 50, 0, 50);
166 if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
167 fProfFracPtMatched = new TProfile("fProfFracPtMatched", "recovered target p_{T} / source p_{T}", 15, -50, 250);
168 fOutputList->Add(fProfFracPtMatched);
169 fProfFracNoMatched = new TProfile("fProfFracNoMatched", "recovered target constituents / source constituents", 15, -50, 250);
170 fOutputList->Add(fProfFracNoMatched);
171 }
172 // the analysis summary histo which stores all the analysis flags is always written to file
173 fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
174 fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
175 fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
176 fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
177 fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
178 fOutputList->Add(fProfQAMatched);
179 fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
180 fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
181 fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
182 fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
183 fOutputList->Add(fProfQA);
184 fProfFracPtJets = new TProfile("fProfFracPtJets", "source p_{T} / target p_{T}", 15, -50, 250);
185 fOutputList->Add(fProfFracPtJets);
186 fProfFracNoJets = new TProfile("fProfFracNoJets", "source constituents / target constituents", 15, -50, 250);
187 fOutputList->Add(fProfFracNoJets);
188 fOutputList->Sort();
189 PostData(1, fOutputList);
190}
191//_____________________________________________________________________________
192TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
193{
194 // book a TH1F and connect it to the output container
195 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
196 if(!fOutputList) return 0x0;
197 TString title(name);
198 title += Form(";%s;[counts]", x);
199 TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
200 histogram->Sumw2();
201 if(append) fOutputList->Add(histogram);
202 return histogram;
203}
204//_____________________________________________________________________________
205TH2F* 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)
206{
207 // book a TH2F and connect it to the output container
208 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
209 if(!fOutputList) return 0x0;
210 TString title(name);
211 title += Form(";%s;%s", x, y);
212 TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
213 histogram->Sumw2();
214 if(append) fOutputList->Add(histogram);
215 return histogram;
216}
217//_____________________________________________________________________________
218Bool_t AliAnalysisTaskJetMatching::Notify()
219{
220 // for each file get the number of trials and pythia cross section
221 // see $ALICE_ROOT/PWGJE/AliAnalysisTaskJetProperties.cxx
222 // $ALICE_ROOT/PWG/Tools/AliAnalysisHelperJetTasks.cxx
223 // this function is implenented here temporarily to avoid introducing a dependency
224 // later on this could just be a call to a static helper function
225 if(!fDoDetectorResponse) return kTRUE;
226 Float_t xsection(0), ftrials(1);
227 fAvgTrials = ftrials;
228 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
229 if(tree) {
230 TFile *curfile = tree->GetCurrentFile();
231 if (!curfile) return kFALSE;
232 TString file(curfile->GetName());
233 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),"");
234 else file.ReplaceAll(gSystem->BaseName(file.Data()),"");
235 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
236 if(!fxsec) {
237 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
238 if(fxsec) {
239 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
240 if(key) {
241 TList *list = dynamic_cast<TList*>(key->ReadObj());
242 if(list) {
243 xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
244 ftrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
245 }
246 }
247 fxsec->Close();
248 }
249 } else {
250 TTree *xtree = (TTree*)fxsec->Get("Xsection");
251 if(xtree) {
252 UInt_t _ftrials = 0;
253 Double_t _xsection = 0;
254 xtree->SetBranchAddress("xsection",&_xsection);
255 xtree->SetBranchAddress("ntrials",&_ftrials);
256 xtree->GetEntry(0);
257 ftrials = _ftrials;
258 xsection = _xsection;
259 }
260 fxsec->Close();
261 }
262 fh1Xsec->Fill("<#sigma>",xsection);
263 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
264 if(ftrials >= nEntries && nEntries > 0.) fAvgTrials = ftrials/nEntries;
265 fh1Trials->Fill("#sum{ntrials}",ftrials);
266 }
267 return kTRUE;
268}
269//_____________________________________________________________________________
270Bool_t AliAnalysisTaskJetMatching::Run()
271{
272 // execute once for each event
273 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
274 if(!(InputEvent() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
275 if(fh1AvgTrials) fh1AvgTrials->Fill("#sum{avg ntrials}", fAvgTrials);
276 // step one: do geometric matching
277 switch (fMatchingScheme) {
278 case kGeoEtaPhi : {
279 DoGeometricMatchingEtaPhi();
280 } break;
281 case kGeoR : {
282 DoGeometricMatchingR();
283 } break;
284 default : break;
285 }
286 // break if no jet was matched
287 if(!fMatchedJetContainer[1][0]) return kTRUE;
288 // optional step two: get a bijection (avoid duplicate matches)
289 if(fGetBijection) GetBijection();
290 // optional step three: match constituents within matched jets
291 if(fMatchConstituents) DoConstituentMatching();
292 // stream data to output
293 PostMatchedJets();
294 FillMatchedJetHistograms();
295 if(fDebug > 0) PrintInfo();
296 PostData(1, fOutputList);
297 return kTRUE;
298}
299//_____________________________________________________________________________
300void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
301{
302 // do geometric matching based on eta phi distance between jet centers
303 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
304 fNoMatchedJets = 0; // reset the matched jet counter
305 Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
306 for(Int_t i(0); i < iSource; i++) {
307 AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
308 if(!PassesCuts(sourceJet, 0)) continue;
309 for(Int_t j(0); j < iTarget; j++) {
310 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
311 if(!PassesCuts(targetJet, 1)) continue;
312 if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
313 if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
314 Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
315 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
316 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
317 if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) { // accept the jets as matching
318 Bool_t isBestMatch(kTRUE);
319 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
320 if(fDebug > 0) printf(" > Entering first bbijection test \n");
321 for(Int_t k(i); k < iSource; k++) {
322 AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
323 if(PassesCuts(candidateSourceJet, 0)) continue;
324 if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
325 if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
326 isBestMatch = kFALSE;
327 break;
328 }
329 }
330 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
331 }
332 if(isBestMatch) {
333 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
334 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
335 fNoMatchedJets++;
336 }
337 if(fNoMatchedJets > 199) { // maximum to keep the cache at reasonable size
338 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
339 return;
340 }
341 }
342 }
343 }
344 }
345}
346//_____________________________________________________________________________
347void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
348{
349 // do geometric matching based on shortest path between jet centers
350 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
351 fNoMatchedJets = 0; // reset the matched jet counter
352 Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
353 for(Int_t i(0); i < iSource; i++) {
354 AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
355 if(!PassesCuts(sourceJet, 0)) continue;
356 for(Int_t j(0); j < iTarget; j++) {
357 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
358 if(!PassesCuts(targetJet, 1)) continue;
359 else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
360 if(GetR(sourceJet, targetJet) <= fMatchR) {
361 Bool_t isBestMatch(kTRUE);
362 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
363 if(fDebug > 0) printf(" > Entering first bijection test \n");
364 for(Int_t k(i); k < iSource; k++) {
365 AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
366 if(!PassesCuts(candidateSourceJet, 0)) continue;
367 if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
368 if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
369 isBestMatch = kFALSE;
370 break;
371 }
372 }
373 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
374 }
375 if(isBestMatch) {
376 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
377 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
378 fNoMatchedJets++;
379 }
380 if(fNoMatchedJets > 99) {
381 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
382 return;
383 }
384 }
385 }
386 }
387}
388//_____________________________________________________________________________
389void AliAnalysisTaskJetMatching::DoConstituentMatching()
390{
391 // match constituents within matched jets
392 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
393 if(!fTracks) {
394 AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
395 return; // coverity ...
396 }
397 for(Int_t i(0); i < fNoMatchedJets; i++) {
398 AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
399 AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
400 if(sourceJet && targetJet) { // duplicate check: slot migth be NULL
401 Double_t targetPt(0);
402 Int_t iSJ(sourceJet->GetNumberOfTracks());
403 Int_t iTJ(targetJet->GetNumberOfTracks());
404 Int_t overlap(0), alreadyFound(0);
405 for(Int_t j(0); j < iSJ; j++) {
406 alreadyFound = 0;
407 Int_t idSource((Int_t)sourceJet->TrackAt(j));
408 for(Int_t k(0); k < iTJ; k++) { // compare all target tracks to the source track
409 if(idSource == targetJet->TrackAt(k) && alreadyFound == 0) {
410 overlap++;
411 alreadyFound++; // avoid possible duplicate matching
412 AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
413 if(vp) targetPt += vp->Pt();
414 continue;
415 }
416 }
417 }
418 if((float)overlap/(float)iSJ < fMinFracRecoveredConstituents || targetPt/sourceJet->Pt() < fMinFracRecoveredConstituentPt) {
419 if(fDebug > 0) printf(" \n > Purging jet, recovered constituents ratio %i / %i = %.2f < or pt ratio %.2f / %.2f = %.2f < %.2f",
420 overlap, iSJ, (float)overlap/(float)iSJ, targetPt, sourceJet->Pt(), targetPt/sourceJet->Pt(), fMinFracRecoveredConstituentPt);
421 fMatchedJetContainer[i][0] = 0x0;
422 fMatchedJetContainer[i][1] = 0x0;
423 continue;
424 }
425 if(sourceJet->Pt() > 0) {
426 Double_t sourceRho(0), targetRho(0);
427 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(sourceJet->Phi(), fSourceRadius)*sourceJet->Area();
428 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(targetJet->Phi(), fTargetRadius)*targetJet->Area();
429 fProfFracPtMatched->Fill(sourceJet->Pt()-sourceRho, (targetPt-targetRho) / (sourceJet->Pt()-sourceRho));
430 fProfFracPtJets->Fill(sourceJet->Pt()-sourceRho, (targetJet->Pt()-targetRho) / (sourceJet->Pt()-sourceRho));
431 fProfFracNoMatched->Fill(sourceJet->Pt()-sourceRho, (double)overlap / (double)sourceJet->GetNumberOfTracks());
432 fProfFracNoJets->Fill(sourceJet->Pt()-sourceRho, (double)targetJet->GetNumberOfTracks() / (double)sourceJet->GetNumberOfTracks());
433 }
434 if(fDebug > 0) {
435 printf("\n > Jet A: %i const\t", iSJ);
436 printf(" > Jet B %i const\t", iTJ);
437 printf(" -> OVERLAP: %i tracks <- \n", overlap);
438 }
439 }
440 }
441}
442//_____________________________________________________________________________
443void AliAnalysisTaskJetMatching::GetBijection()
444{
445 // bijection of source and matched jets, based on closest distance between jets
446 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
447 for(Int_t i(0); i < fNoMatchedJets; i++) {
448 for(Int_t j(i+1); j < fNoMatchedJets; j++) {
449 if(fMatchedJetContainer[i][0] == fMatchedJetContainer[j][0]) {
450 // found source with two targets, now see which target is closer to the source
451 if(!(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1] && fMatchedJetContainer[j][0] && fMatchedJetContainer[j][1] )) continue;
452 Double_t rA(GetR(fMatchedJetContainer[i][0], fMatchedJetContainer[i][1])); // distance between connection A = i
453 Double_t rB(GetR(fMatchedJetContainer[j][0], fMatchedJetContainer[j][1])); // distance between connection B = j
454 if (rB > rA) { // jet two is far away, clear it from both target and source list
455 fMatchedJetContainer[j][0] = 0x0;
456 fMatchedJetContainer[j][1] = 0x0;
457 } else { // jet one is far away, clear it from both target and source list
458 fMatchedJetContainer[i][0] = fMatchedJetContainer[j][0]; // switch to avoid breaking loop
459 fMatchedJetContainer[i][1] = fMatchedJetContainer[j][1];
460 fMatchedJetContainer[j][0] = 0x0; // clear duplicate jet from cache
461 fMatchedJetContainer[j][1] = 0x0;
462 }
463 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
464 }
465 }
466 }
467}
468//_____________________________________________________________________________
469void AliAnalysisTaskJetMatching::PostMatchedJets()
470{
471 // post matched jets
472 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
473 fMatchedJets->Delete(); // should be a NULL operation, but added just in case
474 for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
475 if(fMatchedJetContainer[i][1]) { // duplicate jets will have NULL value here and are skipped
476 new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
477 p++;
478 }
479 }
480}
481//_____________________________________________________________________________
482void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
483{
484 // fill the analysis summary histrogram, saves all relevant analysis settigns
485 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
486 fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
487 fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
488 fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
489 fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
490 fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
491 fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
492 fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
493 fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
494 fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
495 fHistAnalysisSummary->SetBinContent(5, fMatchEta);
496 fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
497 fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
498 fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
499 fHistAnalysisSummary->SetBinContent(7, fMatchR);
500 fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
501 fHistAnalysisSummary->SetBinContent(8, 1.);
502}
503//_____________________________________________________________________________
504void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
505{
506 // fill matched jet histograms
507 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
508 for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
509 AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
510 if(!source) continue;
511 else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
512 Double_t sourceRho(0), targetRho(0);
513 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
514 fHistSourceJetPt->Fill(source->Pt()-sourceRho);
515 fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
516 for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
517 AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
518 if(target) {
519 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
520 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
521 fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));
522 fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
523 fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
524 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
525 if(j==0) {
526 fHistTargetJetPt->Fill(target->Pt()-targetRho);
527 fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
528 }
529 }
530 }
531 }
532 for(Int_t i(0); i < fNoMatchedJets; i++) {
533 if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
534 Double_t sourceRho(0), targetRho(0);
535 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
536 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
537 fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
538 fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
539 fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
540 fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
541 fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
542 fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
543 fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
544
545 fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
546 if(fDoDetectorResponse) {
547 fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
548 fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
549 fHistDetectorResponseProb->Fill((fMatchedJetContainer[i][1]->Pt()-fMatchedJetContainer[i][0]->Pt())/fMatchedJetContainer[i][0]->Pt(), fMatchedJetContainer[i][0]->Pt());
550 }
551 }
552 }
553}
554//_____________________________________________________________________________
555void AliAnalysisTaskJetMatching::PrintInfo() const
556{
557 // print some info
558 printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
559 printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
560 printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
561 if(fDebug > 3) InputEvent()->GetList()->ls();
562}
563//_____________________________________________________________________________
564void AliAnalysisTaskJetMatching::Terminate(Option_t *)
565{
566 // terminate
567}
568//_____________________________________________________________________________