2 // General task to match two sets of jets
4 // This task takes two TClonesArray's as input:
5 // [1] fSourceJets - e.g. pythia jets
6 // [2] fTargetJets - e.g. a samle containing pythia jets embedded in a pbpb event
7 // the task will try to match jets from the source array to the target array, and
8 // save the found TARGET jets in a new array 'fMatchedJets'
9 // the purpose of this task is twofold
10 // [1] matching for embedding studies
11 // [2] matching to create a detector response function
13 // matching is done in three steps
14 // [1] geometric matching, where jets are matched by R = sqrt(dphi*dphi-deta*deta) or directly via dphi and deta
15 // [2] optional injection / bijection check
16 // in this check, it is made sure that fSourceJets -> fMatchedJets is either an injective non-surjective
17 // or bijective map, depending on the matching resolution which is chosen for step [1]
18 // so that each source jet has a unique target and vice-versa.
19 // if R (or dphi, deta) are proportional to, or larger than, the jet radius, matching tends to be
20 // bijective (each source has a target), if R is chosen to be very small, source jets might be lost
21 // as no target can be found.
22 // the mapping is done in such a way that each target is matched to its closest source and each source
23 // is mapped to its closest target
24 // [3] constituent matching
25 // - how many constituents of the source jet are present in the matched jet?
26 // a cut on the constituent fraction can be performed (not recommended)
27 // - how much of the original pt is recovered in the matched jet?
28 // a cut on the fraction of recovered / original pt can be performed (recommended)
31 // to obtain a detector respose function, supply
32 // [1] fSourceJets - particle level jets
33 // [2] fTargetJets - detector level jets
35 // Author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
36 // rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
39 #include <TClonesArray.h>
48 #include <AliAnalysisTask.h>
49 #include <AliAnalysisManager.h>
51 #include <AliVEvent.h>
52 #include <AliVParticle.h>
53 // emcal jet framework includes
54 #include <AliEmcalJet.h>
55 #include <AliAnalysisTaskJetMatching.h>
56 #include <AliLocalRhoParameter.h>
58 class AliAnalysisTaskJetMatching;
61 ClassImp(AliAnalysisTaskJetMatching)
63 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", kTRUE),
64 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), 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) {
65 // default constructor
66 ClearMatchedJetsCache();
68 //_____________________________________________________________________________
69 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
70 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), 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) {
72 ClearMatchedJetsCache();
73 DefineInput(0, TChain::Class());
74 DefineOutput(1, TList::Class());
76 //_____________________________________________________________________________
77 AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
80 if(fOutputList) delete fOutputList;
82 //_____________________________________________________________________________
83 void AliAnalysisTaskJetMatching::ExecOnce()
85 // initialize the anaysis
86 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
87 // get the stand alone jets from the input event
88 fSourceJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSourceJetsName.Data()));
89 if(!fSourceJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
90 fTargetJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTargetJetsName.Data()));
91 if(!fTargetJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
92 // append the list of matched jets to the event
93 fMatchedJets->Delete();
94 if (!(InputEvent()->FindListObject(fMatchedJetsName))) InputEvent()->AddObject(fMatchedJets);
95 else AliFatal(Form("%s: Object with name %s already in event! Aborting", GetName(), fMatchedJetsName.Data()));
96 FillAnalysisSummaryHistogram();
98 case kSourceLocalRho : {
99 fSourceRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fSourceRhoName));
100 if(!fSourceRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fSourceRhoName.Data()));
104 switch (fTargetBKG) {
105 case kTargetLocalRho : {
106 fTargetRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fTargetRhoName));
107 if(!fTargetRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fTargetRhoName.Data()));
111 AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
112 if(fDoDetectorResponse) fMatchConstituents = kFALSE;
114 //_____________________________________________________________________________
115 void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
117 // create output objects
118 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
119 // add the matched jets array to the event
120 fMatchedJets = new TClonesArray("AliEmcalJet");
121 fMatchedJets->SetName(fMatchedJetsName);
122 fOutputList = new TList();
123 fOutputList->SetOwner(kTRUE);
124 // add analysis histograms
125 fHistUnsortedCorrelation = BookTH1F("fHistUnsortedCorrelation", "#Delta #varphi unsorted", 50, 0, TMath::Pi());
126 fHistMatchedCorrelation = BookTH1F("fHistMatchedCorrelation", "#Delta #varphi matched", 50, 0, TMath::Pi());
127 fHistSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
128 fHistMatchedSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistMatchedParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistMatchedSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
129 fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
130 fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
131 fHistSourceMatchedJetPt = (fDoDetectorResponse) ? BookTH2F("fHistPartDetJetPt", "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);
132 fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents", 50, 0, 50);
133 fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents", 50, 0, 50);
134 fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents", 50, 0, 50);
135 if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
136 fProfFracPtMatched = new TProfile("fProfFracPtMatched", "fProfFracPtMatched", 15, -50, 250);
137 fOutputList->Add(fProfFracPtMatched);
138 fProfFracNoMatched = new TProfile("fProfFracNoMatched", "fProfFracNoMatched", 15, -50, 250);
139 fOutputList->Add(fProfFracNoMatched);
141 // the analysis summary histo which stores all the analysis flags is always written to file
142 fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
143 fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
144 fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
145 fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
146 fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
147 fOutputList->Add(fProfQAMatched);
148 fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
149 fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
150 fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
151 fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
152 fOutputList->Add(fProfQA);
153 fProfFracPtJets = new TProfile("fProfFracPtJets", "fProfFracPtJets", 15, -50, 250);
154 fOutputList->Add(fProfFracPtJets);
155 fProfFracNoJets = new TProfile("fProfFracNoJets", "fProfFracNoJets", 15, -50, 250);
156 fOutputList->Add(fProfFracNoJets);
158 PostData(1, fOutputList);
160 //_____________________________________________________________________________
161 TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
163 // book a TH1F and connect it to the output container
164 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
165 if(!fOutputList) return 0x0;
167 title += Form(";%s;[counts]", x);
168 TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
170 if(append) fOutputList->Add(histogram);
173 //_____________________________________________________________________________
174 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)
176 // book a TH2F and connect it to the output container
177 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
178 if(!fOutputList) return 0x0;
180 title += Form(";%s;%s", x, y);
181 TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
183 if(append) fOutputList->Add(histogram);
186 //_____________________________________________________________________________
187 Bool_t AliAnalysisTaskJetMatching::Run()
189 // execute once for each event
190 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
191 if(!(InputEvent() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
192 // step one: do geometric matching
193 switch (fMatchingScheme) {
195 DoGeometricMatchingEtaPhi();
198 DoGeometricMatchingR();
202 // break if no jet was matched
203 if(!fMatchedJetContainer[1][0]) return kTRUE;
204 // optional step two: get a bijection (avoid duplicate matches)
205 if(fGetBijection) GetBijection();
206 // optional step three: match constituents within matched jets
207 if(fMatchConstituents) DoConstituentMatching();
208 // stream data to output
210 FillMatchedJetHistograms();
211 if(fDebug > 0) PrintInfo();
212 PostData(1, fOutputList);
215 //_____________________________________________________________________________
216 void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
218 // do geometric matching based on eta phi distance between jet centers
219 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
220 fNoMatchedJets = 0; // reset the matched jet counter
221 Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
222 for(Int_t i(0); i < iSource; i++) {
223 AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
224 if(!PassesCuts(sourceJet)) continue;
225 if(fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
226 for(Int_t j(0); j < iTarget; j++) {
227 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
228 if(!PassesCuts(targetJet)) continue;
229 if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
230 if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
231 Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
232 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
233 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
234 if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) { // accept the jets as matching
235 Bool_t isBestMatch(kTRUE);
236 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
237 if(fDebug > 0) printf(" > Entering first bbijection test \n");
238 for(Int_t k(i); k < iSource; k++) {
239 AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
240 if(!PassesCuts(candidateSourceJet)) continue;
241 if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
242 if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
243 if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
244 isBestMatch = kFALSE;
248 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
251 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
252 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
255 if(fNoMatchedJets > 199) { // maximum to keep the cache at reasonable size
256 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
264 //_____________________________________________________________________________
265 void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
267 // do geometric matching based on shortest path between jet centers
268 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
269 fNoMatchedJets = 0; // reset the matched jet counter
270 Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
271 for(Int_t i(0); i < iSource; i++) {
272 AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
273 if(!PassesCuts(sourceJet)) continue;
274 else if (fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
275 for(Int_t j(0); j < iTarget; j++) {
276 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
277 if(!PassesCuts(targetJet)) continue;
278 else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
279 if(GetR(sourceJet, targetJet) <= fMatchR) {
280 Bool_t isBestMatch(kTRUE);
281 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
282 if(fDebug > 0) printf(" > Entering first bijection test \n");
283 for(Int_t k(i); k < iSource; k++) {
284 AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
285 if(!PassesCuts(candidateSourceJet)) continue;
286 if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
287 if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
288 if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
289 isBestMatch = kFALSE;
293 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
296 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
297 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
300 if(fNoMatchedJets > 99) {
301 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
308 //_____________________________________________________________________________
309 void AliAnalysisTaskJetMatching::DoConstituentMatching()
311 // match constituents within matched jets
312 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
314 AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
315 return; // coverity ...
317 for(Int_t i(0); i < fNoMatchedJets; i++) {
318 AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
319 AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
320 if(sourceJet && targetJet) { // duplicate check: slot migth be NULL
321 Double_t targetPt(0);
322 Int_t iSJ(sourceJet->GetNumberOfTracks());
323 Int_t iTJ(targetJet->GetNumberOfTracks());
324 Int_t overlap(0), alreadyFound(0);
325 for(Int_t j(0); j < iSJ; j++) {
327 Int_t idSource((Int_t)sourceJet->TrackAt(j));
328 for(Int_t k(0); k < iTJ; k++) { // compare all target tracks to the source track
329 if(idSource == targetJet->TrackAt(k) && alreadyFound == 0) {
331 alreadyFound++; // avoid possible duplicate matching
332 AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
333 if(vp) targetPt += vp->Pt();
338 if((float)overlap/(float)iSJ < fMinFracRecoveredConstituents || targetPt/sourceJet->Pt() < fMinFracRecoveredConstituentPt) {
339 if(fDebug > 0) printf(" \n > Purging jet, recovered constituents ratio %i / %i = %.2f < or pt ratio %.2f / %.2f = %.2f < %.2f",
340 overlap, iSJ, (float)overlap/(float)iSJ, targetPt, sourceJet->Pt(), targetPt/sourceJet->Pt(), fMinFracRecoveredConstituentPt);
341 fMatchedJetContainer[i][0] = 0x0;
342 fMatchedJetContainer[i][1] = 0x0;
345 if(sourceJet->Pt() > 0) {
346 Double_t sourceRho(0), targetRho(0);
347 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(sourceJet->Phi(), fSourceRadius)*sourceJet->Area();
348 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(targetJet->Phi(), fTargetRadius)*targetJet->Area();
349 fProfFracPtMatched->Fill(sourceJet->Pt()-sourceRho, (targetPt-targetRho) / (sourceJet->Pt()-sourceRho));
350 fProfFracPtJets->Fill(sourceJet->Pt()-sourceRho, (targetJet->Pt()-targetRho) / (sourceJet->Pt()-sourceRho));
351 fProfFracNoMatched->Fill(sourceJet->Pt()-sourceRho, (double)overlap / (double)sourceJet->GetNumberOfTracks());
352 fProfFracNoJets->Fill(sourceJet->Pt()-sourceRho, (double)targetJet->GetNumberOfTracks() / (double)sourceJet->GetNumberOfTracks());
355 printf("\n > Jet A: %i const\t", iSJ);
356 printf(" > Jet B %i const\t", iTJ);
357 printf(" -> OVERLAP: %i tracks <- \n", overlap);
362 //_____________________________________________________________________________
363 void AliAnalysisTaskJetMatching::GetBijection()
365 // bijection of source and matched jets, based on closest distance between jets
366 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
367 for(Int_t i(0); i < fNoMatchedJets; i++) {
368 for(Int_t j(i+1); j < fNoMatchedJets; j++) {
369 if(fMatchedJetContainer[i][0] == fMatchedJetContainer[j][0]) {
370 // found source with two targets, now see which target is closer to the source
371 if(!(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1] && fMatchedJetContainer[j][0] && fMatchedJetContainer[j][1] )) continue;
372 Double_t rA(GetR(fMatchedJetContainer[i][0], fMatchedJetContainer[i][1])); // distance between connection A = i
373 Double_t rB(GetR(fMatchedJetContainer[j][0], fMatchedJetContainer[j][1])); // distance between connection B = j
374 if (rB > rA) { // jet two is far away, clear it from both target and source list
375 fMatchedJetContainer[j][0] = 0x0;
376 fMatchedJetContainer[j][1] = 0x0;
377 } else { // jet one is far away, clear it from both target and source list
378 fMatchedJetContainer[i][0] = fMatchedJetContainer[j][0]; // switch to avoid breaking loop
379 fMatchedJetContainer[i][1] = fMatchedJetContainer[j][1];
380 fMatchedJetContainer[j][0] = 0x0; // clear duplicate jet from cache
381 fMatchedJetContainer[j][1] = 0x0;
383 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
388 //_____________________________________________________________________________
389 void AliAnalysisTaskJetMatching::PostMatchedJets()
392 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
393 fMatchedJets->Delete(); // should be a NULL operation, but added just in case
394 for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
395 if(fMatchedJetContainer[i][1]) { // duplicate jets will have NULL value here and are skipped
396 new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
401 //_____________________________________________________________________________
402 void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
404 // fill the analysis summary histrogram, saves all relevant analysis settigns
405 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
406 fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
407 fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
408 fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
409 fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
410 fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
411 fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
412 fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
413 fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
414 fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
415 fHistAnalysisSummary->SetBinContent(5, fMatchEta);
416 fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
417 fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
418 fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
419 fHistAnalysisSummary->SetBinContent(7, fMatchR);
420 fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
421 fHistAnalysisSummary->SetBinContent(8, 1.);
423 //_____________________________________________________________________________
424 void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
426 // fill matched jet histograms
427 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
428 for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
429 AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
430 if(!source) continue;
431 else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
432 Double_t sourceRho(0), targetRho(0);
433 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
434 fHistSourceJetPt->Fill(source->Pt()-sourceRho);
435 fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
436 for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
437 AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
439 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
440 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
441 fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));
442 fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
443 fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
444 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
446 fHistTargetJetPt->Fill(target->Pt()-targetRho);
447 fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
452 for(Int_t i(0); i < fNoMatchedJets; i++) {
453 if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
454 Double_t sourceRho(0), targetRho(0);
455 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
456 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
457 fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
458 fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
459 fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
460 fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
461 fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
462 fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
463 fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
464 fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
465 if(fDoDetectorResponse) {
466 fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
467 fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
472 //_____________________________________________________________________________
473 void AliAnalysisTaskJetMatching::PrintInfo() const
476 printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
477 printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
478 printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
479 if(fDebug > 3) InputEvent()->GetList()->ls();
481 //_____________________________________________________________________________
482 void AliAnalysisTaskJetMatching::Terminate(Option_t *)
486 //_____________________________________________________________________________