]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add code for full-charged response (from Filip)
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetMatching.cxx
CommitLineData
f1354962 1//
2// General task to match two sets of jets
3//
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'
c2460f8a 9// the purpose of this task is twofold
10// [1] matching for embedding studies
11// [2] matching to create a detector response function
f1354962 12//
285db795 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)
29//
c2460f8a 30// detector response
31// to obtain a detector respose function, supply
32// [1] fSourceJets - particle level jets
33// [2] fTargetJets - detector level jets
34//
f1354962 35// Author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
36// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
37
38// root includes
39#include <TClonesArray.h>
40#include <TChain.h>
41#include <TMath.h>
42#include <TF1.h>
43#include <TF2.h>
44#include <TH1F.h>
45#include <TH2F.h>
46#include <TProfile.h>
47// aliroot includes
48#include <AliAnalysisTask.h>
49#include <AliAnalysisManager.h>
50#include <AliLog.h>
51#include <AliVEvent.h>
52#include <AliVParticle.h>
53// emcal jet framework includes
54#include <AliEmcalJet.h>
55#include <AliAnalysisTaskJetMatching.h>
f6d1b1a7 56#include <AliLocalRhoParameter.h>
f1354962 57
58class AliAnalysisTaskJetMatching;
59using namespace std;
60
61ClassImp(AliAnalysisTaskJetMatching)
62
e9c2d4f3 63AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJetDev("AliAnalysisTaskJetMatching", kTRUE),
2e3ee924 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) {
f1354962 65 // default constructor
66 ClearMatchedJetsCache();
67}
68//_____________________________________________________________________________
e9c2d4f3 69AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJetDev(name, kTRUE),
2e3ee924 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) {
f1354962 71 // constructor
72 ClearMatchedJetsCache();
73 DefineInput(0, TChain::Class());
74 DefineOutput(1, TList::Class());
75}
76//_____________________________________________________________________________
77AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
78{
79 // destructor
80 if(fOutputList) delete fOutputList;
81}
82//_____________________________________________________________________________
83void AliAnalysisTaskJetMatching::ExecOnce()
84{
85 // initialize the anaysis
86 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
f1354962 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();
f6d1b1a7 97 switch (fSourceBKG) {
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()));
101 } break;
102 default : break;
103 }
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()));
108 } break;
109 default : break;
110 }
e9c2d4f3 111 AliAnalysisTaskEmcalJetDev::ExecOnce(); // init base class
c2460f8a 112 if(fDoDetectorResponse) fMatchConstituents = kFALSE;
f1354962 113}
114//_____________________________________________________________________________
115void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
116{
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());
2e3ee924 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);
c2460f8a 129 fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
2e3ee924 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]", 150, -50, 250, 150, -50, 250) : BookTH2F("fHistSourceMatchedJetPt", "source jet p_{t} [GeV/c]", "matched jet p_{t} [GeV/c]", 150, -50, 250, 150, -50, 250);
f1354962 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);
2e3ee924 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);
140 }
f1354962 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);
c2460f8a 153 fProfFracPtJets = new TProfile("fProfFracPtJets", "fProfFracPtJets", 15, -50, 250);
0aaf31ab 154 fOutputList->Add(fProfFracPtJets);
c2460f8a 155 fProfFracNoJets = new TProfile("fProfFracNoJets", "fProfFracNoJets", 15, -50, 250);
0aaf31ab 156 fOutputList->Add(fProfFracNoJets);
f1354962 157 fOutputList->Sort();
158 PostData(1, fOutputList);
159}
160//_____________________________________________________________________________
161TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
162{
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;
166 TString title(name);
167 title += Form(";%s;[counts]", x);
168 TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
169 histogram->Sumw2();
170 if(append) fOutputList->Add(histogram);
171 return histogram;
172}
173//_____________________________________________________________________________
174TH2F* 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)
175{
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;
179 TString title(name);
180 title += Form(";%s;%s", x, y);
181 TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
182 histogram->Sumw2();
183 if(append) fOutputList->Add(histogram);
184 return histogram;
185}
186//_____________________________________________________________________________
187Bool_t AliAnalysisTaskJetMatching::Run()
188{
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() && fSourceJetsName && fTargetJets && IsEventSelected())) return kFALSE;
285db795 192 // step one: do geometric matching
f1354962 193 switch (fMatchingScheme) {
194 case kGeoEtaPhi : {
195 DoGeometricMatchingEtaPhi();
196 } break;
197 case kGeoR : {
198 DoGeometricMatchingR();
199 } break;
285db795 200 default : break;
f1354962 201 }
285db795 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();
f1354962 208 // stream data to output
209 PostMatchedJets();
210 FillMatchedJetHistograms();
211 if(fDebug > 0) PrintInfo();
212 PostData(1, fOutputList);
213 return kTRUE;
214}
215//_____________________________________________________________________________
285db795 216void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
f1354962 217{
285db795 218 // do geometric matching based on eta phi distance between jet centers
f1354962 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;
f6d1b1a7 225 if(fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
f1354962 226 for(Int_t j(0); j < iTarget; j++) {
227 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
228 if(!PassesCuts(targetJet)) continue;
f6d1b1a7 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();
285db795 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;
245 break;
246 }
247 }
248 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
249 }
250 if(isBestMatch) {
251 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
252 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
253 fNoMatchedJets++;
254 }
c2460f8a 255 if(fNoMatchedJets > 199) { // maximum to keep the cache at reasonable size
f6d1b1a7 256 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
257 return;
258 }
f1354962 259 }
260 }
261 }
262 }
263}
264//_____________________________________________________________________________
285db795 265void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
f1354962 266{
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;
f6d1b1a7 274 else if (fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
f1354962 275 for(Int_t j(0); j < iTarget; j++) {
276 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
277 if(!PassesCuts(targetJet)) continue;
f6d1b1a7 278 else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
279 if(GetR(sourceJet, targetJet) <= fMatchR) {
285db795 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;
290 break;
291 }
292 }
293 if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
294 }
295 if(isBestMatch) {
296 fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
297 fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
298 fNoMatchedJets++;
299 }
f1354962 300 if(fNoMatchedJets > 99) {
301 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
302 return;
303 }
304 }
305 }
306 }
307}
308//_____________________________________________________________________________
285db795 309void AliAnalysisTaskJetMatching::DoConstituentMatching()
f1354962 310{
285db795 311 // match constituents within matched jets
f1354962 312 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
313 if(!fTracks) {
314 AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
315 return; // coverity ...
316 }
317 for(Int_t i(0); i < fNoMatchedJets; i++) {
0aaf31ab 318 AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
319 AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
f1354962 320 if(sourceJet && targetJet) { // duplicate check: slot migth be NULL
0aaf31ab 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++) {
326 alreadyFound = 0;
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) {
330 overlap++;
331 alreadyFound++; // avoid possible duplicate matching
332 AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
333 if(vp) targetPt += vp->Pt();
334 continue;
335 }
f1354962 336 }
337 }
285db795 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);
f6d1b1a7 341 fMatchedJetContainer[i][0] = 0x0;
342 fMatchedJetContainer[i][1] = 0x0;
343 continue;
344 }
0aaf31ab 345 if(sourceJet->Pt() > 0) {
f6d1b1a7 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());
0aaf31ab 353 }
354 if(fDebug > 0) {
285db795 355 printf("\n > Jet A: %i const\t", iSJ);
356 printf(" > Jet B %i const\t", iTJ);
357 printf(" -> OVERLAP: %i tracks <- \n", overlap);
0aaf31ab 358 }
f1354962 359 }
360 }
361}
362//_____________________________________________________________________________
285db795 363void AliAnalysisTaskJetMatching::GetBijection()
f1354962 364{
285db795 365 // bijection of source and matched jets, based on closest distance between jets
f1354962 366 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
f1354962 367 for(Int_t i(0); i < fNoMatchedJets; i++) {
368 for(Int_t j(i+1); j < fNoMatchedJets; j++) {
f6d1b1a7 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;
285db795 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
f6d1b1a7 375 fMatchedJetContainer[j][0] = 0x0;
376 fMatchedJetContainer[j][1] = 0x0;
285db795 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
f6d1b1a7 381 fMatchedJetContainer[j][1] = 0x0;
f1354962 382 }
f6d1b1a7 383 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
f1354962 384 }
385 }
386 }
387}
388//_____________________________________________________________________________
389void AliAnalysisTaskJetMatching::PostMatchedJets()
390{
391 // post matched jets
392 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
393 for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
394 if(fMatchedJetContainer[i][1]) { // duplicate jets will have NULL value here and are skipped
395 new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
396 p++;
397 }
398 }
399}
400//_____________________________________________________________________________
401void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
402{
403 // fill the analysis summary histrogram, saves all relevant analysis settigns
404 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
405 fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
406 fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
407 fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
408 fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
409 fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
410 fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
411 fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
412 fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
285db795 413 fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
414 fHistAnalysisSummary->SetBinContent(5, fMatchEta);
415 fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
416 fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
417 fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
418 fHistAnalysisSummary->SetBinContent(7, fMatchR);
419 fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
420 fHistAnalysisSummary->SetBinContent(8, 1.);
f1354962 421}
422//_____________________________________________________________________________
c2460f8a 423void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
f1354962 424{
425 // fill matched jet histograms
426 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
427 for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
428 AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
429 if(!source) continue;
c2460f8a 430 else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
f6d1b1a7 431 Double_t sourceRho(0), targetRho(0);
432 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
433 fHistSourceJetPt->Fill(source->Pt()-sourceRho);
f1354962 434 fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
435 for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
436 AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
437 if(target) {
c2460f8a 438 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
f6d1b1a7 439 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
440 fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));
441 fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
442 fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
f1354962 443 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
444 if(j==0) {
f6d1b1a7 445 fHistTargetJetPt->Fill(target->Pt()-targetRho);
f1354962 446 fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
447 }
448 }
449 }
450 }
451 for(Int_t i(0); i < fNoMatchedJets; i++) {
452 if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
f6d1b1a7 453 Double_t sourceRho(0), targetRho(0);
454 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
968e3e0d 455 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
f1354962 456 fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
2e3ee924 457 fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
f6d1b1a7 458 fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
459 fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
460 fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
f1354962 461 fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
462 fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
c2460f8a 463 fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
2e3ee924 464 if(fDoDetectorResponse) {
465 fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
466 fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
467 }
f1354962 468 }
469 }
470}
471//_____________________________________________________________________________
472void AliAnalysisTaskJetMatching::PrintInfo() const
473{
474 // print some info
285db795 475 printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
f1354962 476 printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
477 printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
478 if(fDebug > 3) InputEvent()->GetList()->ls();
479}
480//_____________________________________________________________________________
481void AliAnalysisTaskJetMatching::Terminate(Option_t *)
482{
483 // terminate
484}
485//_____________________________________________________________________________