]> 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
Removing printout (Davide)
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetMatching.cxx
CommitLineData
7f7120c5 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
f1354962 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'
c2460f8a 24// the purpose of this task is twofold
25// [1] matching for embedding studies
26// [2] matching to create a detector response function
f1354962 27//
285db795 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//
c2460f8a 45// detector response
46// to obtain a detector respose function, supply
47// [1] fSourceJets - particle level jets
48// [2] fTargetJets - detector level jets
49//
f1354962 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>
7ec64e26 62#include <TFile.h>
63#include <TTree.h>
64#include <TKey.h>
65#include <TSystem.h>
f1354962 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>
f6d1b1a7 75#include <AliLocalRhoParameter.h>
f1354962 76
77class AliAnalysisTaskJetMatching;
78using namespace std;
79
80ClassImp(AliAnalysisTaskJetMatching)
81
9239b066 82AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", kTRUE),
7ec64e26 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) {
f1354962 84 // default constructor
85 ClearMatchedJetsCache();
86}
87//_____________________________________________________________________________
9239b066 88AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
7ec64e26 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) {
f1354962 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__);
f1354962 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();
f6d1b1a7 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 }
9239b066 130 AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
c2460f8a 131 if(fDoDetectorResponse) fMatchConstituents = kFALSE;
f1354962 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());
2e3ee924 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);
c2460f8a 148 fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
2e3ee924 149 fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
7f7120c5 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);
cd5876d1 151 if(fDoDetectorResponse) {
7ec64e26 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);
cd5876d1 162 }
7f7120c5 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);
2e3ee924 166 if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
7f7120c5 167 fProfFracPtMatched = new TProfile("fProfFracPtMatched", "recovered target p_{T} / source p_{T}", 15, -50, 250);
2e3ee924 168 fOutputList->Add(fProfFracPtMatched);
7f7120c5 169 fProfFracNoMatched = new TProfile("fProfFracNoMatched", "recovered target constituents / source constituents", 15, -50, 250);
2e3ee924 170 fOutputList->Add(fProfFracNoMatched);
171 }
f1354962 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);
7f7120c5 184 fProfFracPtJets = new TProfile("fProfFracPtJets", "source p_{T} / target p_{T}", 15, -50, 250);
0aaf31ab 185 fOutputList->Add(fProfFracPtJets);
7f7120c5 186 fProfFracNoJets = new TProfile("fProfFracNoJets", "source constituents / target constituents", 15, -50, 250);
0aaf31ab 187 fOutputList->Add(fProfFracNoJets);
f1354962 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//_____________________________________________________________________________
7ec64e26 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//_____________________________________________________________________________
f1354962 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__);
459e1638 274 if(!(InputEvent() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
7ec64e26 275 if(fh1AvgTrials) fh1AvgTrials->Fill("#sum{avg ntrials}", fAvgTrials);
285db795 276 // step one: do geometric matching
f1354962 277 switch (fMatchingScheme) {
278 case kGeoEtaPhi : {
279 DoGeometricMatchingEtaPhi();
280 } break;
281 case kGeoR : {
282 DoGeometricMatchingR();
283 } break;
285db795 284 default : break;
f1354962 285 }
285db795 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();
f1354962 292 // stream data to output
293 PostMatchedJets();
294 FillMatchedJetHistograms();
295 if(fDebug > 0) PrintInfo();
296 PostData(1, fOutputList);
297 return kTRUE;
298}
299//_____________________________________________________________________________
285db795 300void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
f1354962 301{
285db795 302 // do geometric matching based on eta phi distance between jet centers
f1354962 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)));
4f6d96b5 308 if(!PassesCuts(sourceJet, 0)) continue;
f1354962 309 for(Int_t j(0); j < iTarget; j++) {
310 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
4f6d96b5 311 if(!PassesCuts(targetJet, 1)) continue;
f6d1b1a7 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();
285db795 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)));
4f6d96b5 323 if(PassesCuts(candidateSourceJet, 0)) continue;
285db795 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 }
c2460f8a 337 if(fNoMatchedJets > 199) { // maximum to keep the cache at reasonable size
f6d1b1a7 338 AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
339 return;
340 }
f1354962 341 }
342 }
343 }
344 }
345}
346//_____________________________________________________________________________
285db795 347void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
f1354962 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)));
4f6d96b5 355 if(!PassesCuts(sourceJet, 0)) continue;
f1354962 356 for(Int_t j(0); j < iTarget; j++) {
357 AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
4f6d96b5 358 if(!PassesCuts(targetJet, 1)) continue;
f6d1b1a7 359 else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
360 if(GetR(sourceJet, targetJet) <= fMatchR) {
285db795 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)));
4f6d96b5 366 if(!PassesCuts(candidateSourceJet, 0)) continue;
285db795 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 }
f1354962 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//_____________________________________________________________________________
285db795 389void AliAnalysisTaskJetMatching::DoConstituentMatching()
f1354962 390{
285db795 391 // match constituents within matched jets
f1354962 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++) {
0aaf31ab 398 AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
399 AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
f1354962 400 if(sourceJet && targetJet) { // duplicate check: slot migth be NULL
0aaf31ab 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 }
f1354962 416 }
417 }
285db795 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);
f6d1b1a7 421 fMatchedJetContainer[i][0] = 0x0;
422 fMatchedJetContainer[i][1] = 0x0;
423 continue;
424 }
0aaf31ab 425 if(sourceJet->Pt() > 0) {
f6d1b1a7 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());
0aaf31ab 433 }
434 if(fDebug > 0) {
285db795 435 printf("\n > Jet A: %i const\t", iSJ);
436 printf(" > Jet B %i const\t", iTJ);
437 printf(" -> OVERLAP: %i tracks <- \n", overlap);
0aaf31ab 438 }
f1354962 439 }
440 }
441}
442//_____________________________________________________________________________
285db795 443void AliAnalysisTaskJetMatching::GetBijection()
f1354962 444{
285db795 445 // bijection of source and matched jets, based on closest distance between jets
f1354962 446 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
f1354962 447 for(Int_t i(0); i < fNoMatchedJets; i++) {
448 for(Int_t j(i+1); j < fNoMatchedJets; j++) {
f6d1b1a7 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;
285db795 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
f6d1b1a7 455 fMatchedJetContainer[j][0] = 0x0;
456 fMatchedJetContainer[j][1] = 0x0;
285db795 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
f6d1b1a7 461 fMatchedJetContainer[j][1] = 0x0;
f1354962 462 }
f6d1b1a7 463 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
f1354962 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__);
459e1638 473 fMatchedJets->Delete(); // should be a NULL operation, but added just in case
f1354962 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);
285db795 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.);
f1354962 502}
503//_____________________________________________________________________________
c2460f8a 504void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
f1354962 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;
c2460f8a 511 else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
f6d1b1a7 512 Double_t sourceRho(0), targetRho(0);
513 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
514 fHistSourceJetPt->Fill(source->Pt()-sourceRho);
f1354962 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) {
c2460f8a 519 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
f6d1b1a7 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()));
f1354962 524 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
525 if(j==0) {
f6d1b1a7 526 fHistTargetJetPt->Fill(target->Pt()-targetRho);
f1354962 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]) {
f6d1b1a7 534 Double_t sourceRho(0), targetRho(0);
535 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
968e3e0d 536 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
f1354962 537 fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
2e3ee924 538 fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
f6d1b1a7 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)));
f1354962 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()));
4f6d96b5 544
c2460f8a 545 fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
2e3ee924 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());
7ec64e26 549 fHistDetectorResponseProb->Fill((fMatchedJetContainer[i][1]->Pt()-fMatchedJetContainer[i][0]->Pt())/fMatchedJetContainer[i][0]->Pt(), fMatchedJetContainer[i][0]->Pt());
2e3ee924 550 }
f1354962 551 }
552 }
553}
554//_____________________________________________________________________________
555void AliAnalysisTaskJetMatching::PrintInfo() const
556{
557 // print some info
285db795 558 printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
f1354962 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//_____________________________________________________________________________