d23ffcb5c108844374b9af32cec4bec4848117d6
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetMatching.cxx
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'
9 // the purpose of this task is twofold
10 // [1] matching for embedding studies
11 // [2] matching to create a detector response function
12 //
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 //
30 // detector response
31 //     to obtain a detector respose function, supply
32 // [1] fSourceJets - particle level jets
33 // [2] fTargetJets - detector level jets
34 // 
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>
56 #include <AliLocalRhoParameter.h>
57
58 class AliAnalysisTaskJetMatching;
59 using namespace std;
60
61 ClassImp(AliAnalysisTaskJetMatching)
62
63 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJetDev("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), 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();
67 }
68 //_____________________________________________________________________________
69 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJetDev(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), 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) {
71     // constructor
72     ClearMatchedJetsCache();
73     DefineInput(0, TChain::Class());
74     DefineOutput(1, TList::Class());
75 }
76 //_____________________________________________________________________________
77 AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
78 {
79     // destructor
80     if(fOutputList)             delete fOutputList;
81 }
82 //_____________________________________________________________________________
83 void AliAnalysisTaskJetMatching::ExecOnce() 
84 {
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();
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     }
111     AliAnalysisTaskEmcalJetDev::ExecOnce(); // init base class
112     if(fDoDetectorResponse) fMatchConstituents = kFALSE;
113 }
114 //_____________________________________________________________________________
115 void 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());
127     fHistSourceJetPt = BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
128     fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
129     fHistMatchedJetPt = BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
130     fHistSourceMatchedJetPt = BookTH2F("fHistSourceMatchedJetPt", "source jet p_{t} [GeV/c]", "matched jet p_{t} [GeV/c]", 150, -50, 250, 150, -50, 250);
131     fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents", 50, 0, 50);
132     fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents", 50, 0, 50);
133     fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents", 50, 0, 50);
134     // the analysis summary histo which stores all the analysis flags is always written to file
135     fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
136     fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
137     fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
138     fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
139     fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
140     fOutputList->Add(fProfQAMatched);
141     fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
142     fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
143     fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
144     fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
145     fOutputList->Add(fProfQA);
146     fProfFracPtMatched = new TProfile("fProfFracPtMatched", "fProfFracPtMatched", 15, -50, 250);
147     fOutputList->Add(fProfFracPtMatched);
148     fProfFracPtJets = new TProfile("fProfFracPtJets", "fProfFracPtJets", 15, -50, 250);
149     fOutputList->Add(fProfFracPtJets);
150     fProfFracNoMatched = new TProfile("fProfFracNoMatched", "fProfFracNoMatched", 15, -50, 250);
151     fOutputList->Add(fProfFracNoMatched);
152     fProfFracNoJets = new TProfile("fProfFracNoJets", "fProfFracNoJets", 15, -50, 250);
153     fOutputList->Add(fProfFracNoJets);
154     fOutputList->Sort();
155     PostData(1, fOutputList);
156 }
157 //_____________________________________________________________________________
158 TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
159 {
160     // book a TH1F and connect it to the output container
161     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
162     if(!fOutputList) return 0x0;
163     TString title(name);
164     title += Form(";%s;[counts]", x);
165     TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
166     histogram->Sumw2();
167     if(append) fOutputList->Add(histogram);
168     return histogram;   
169 }
170 //_____________________________________________________________________________
171 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)
172 {
173     // book a TH2F and connect it to the output container
174     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
175     if(!fOutputList) return 0x0;
176     TString title(name);
177     title += Form(";%s;%s", x, y);
178     TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
179     histogram->Sumw2();
180     if(append) fOutputList->Add(histogram);
181     return histogram;   
182 }
183 //_____________________________________________________________________________
184 Bool_t AliAnalysisTaskJetMatching::Run()
185 {
186     // execute once for each event
187     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
188     if(!(InputEvent() && fSourceJetsName && fTargetJets && IsEventSelected())) return kFALSE;
189     // step one: do geometric matching 
190     switch (fMatchingScheme) {
191         case kGeoEtaPhi : {
192             DoGeometricMatchingEtaPhi();
193         } break;
194         case kGeoR : {
195             DoGeometricMatchingR();
196         } break;
197         default : break;
198     }
199     // break if no jet was matched
200     if(!fMatchedJetContainer[1][0]) return kTRUE;
201     // optional step two: get a bijection (avoid duplicate matches)
202     if(fGetBijection)           GetBijection();
203     // optional step three: match constituents within matched jets
204     if(fMatchConstituents)      DoConstituentMatching();
205     // stream data to output
206     PostMatchedJets();
207     FillMatchedJetHistograms();
208     if(fDebug > 0) PrintInfo();
209     PostData(1, fOutputList);
210     return kTRUE;
211 }
212 //_____________________________________________________________________________
213 void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
214 {
215     // do geometric matching based on eta phi distance between jet centers
216     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
217     fNoMatchedJets = 0; // reset the matched jet counter
218     Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
219     for(Int_t i(0); i < iSource; i++) {
220         AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
221         if(!PassesCuts(sourceJet)) continue;
222         if(fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
223         for(Int_t j(0); j < iTarget; j++) {
224             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
225             if(!PassesCuts(targetJet)) continue;
226             if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
227             if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
228                 Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
229                 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
230                 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
231                 if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {       // accept the jets as matching 
232                 Bool_t isBestMatch(kTRUE);
233                     if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
234                         if(fDebug > 0) printf(" > Entering first bbijection test \n");
235                         for(Int_t k(i); k < iSource; k++) {
236                             AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
237                             if(!PassesCuts(candidateSourceJet)) continue;
238                             if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
239                             if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
240                             if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
241                                 isBestMatch = kFALSE;
242                                 break;
243                             }
244                         }
245                         if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
246                     }
247                     if(isBestMatch) {
248                         fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
249                         fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
250                         fNoMatchedJets++;
251                     }
252                     if(fNoMatchedJets > 199) {   // maximum to keep the cache at reasonable size
253                         AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
254                         return;
255                     }
256                 }
257             }
258         }
259     }
260 }
261 //_____________________________________________________________________________
262 void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
263 {
264     // do geometric matching based on shortest path between jet centers
265     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
266     fNoMatchedJets = 0; // reset the matched jet counter
267     Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
268     for(Int_t i(0); i < iSource; i++) {
269         AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
270         if(!PassesCuts(sourceJet)) continue;
271         else if (fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
272         for(Int_t j(0); j < iTarget; j++) {
273             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
274             if(!PassesCuts(targetJet)) continue;
275             else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
276             if(GetR(sourceJet, targetJet) <= fMatchR) {
277                 Bool_t isBestMatch(kTRUE);
278                 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
279                     if(fDebug > 0) printf(" > Entering first bijection test \n");
280                     for(Int_t k(i); k < iSource; k++) {
281                         AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
282                         if(!PassesCuts(candidateSourceJet)) continue;
283                         if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
284                         if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
285                         if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
286                             isBestMatch = kFALSE;
287                             break;
288                         }
289                     }
290                     if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
291                 }
292                 if(isBestMatch) {
293                     fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
294                     fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
295                     fNoMatchedJets++;
296                 }
297                 if(fNoMatchedJets > 99) {
298                     AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
299                     return;
300                 }
301             }
302         }
303     }
304 }
305 //_____________________________________________________________________________
306 void AliAnalysisTaskJetMatching::DoConstituentMatching()
307 {
308     // match constituents within matched jets 
309     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
310     if(!fTracks) {
311         AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
312         return; // coverity ...
313     }
314     for(Int_t i(0); i < fNoMatchedJets; i++) {
315         AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
316         AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
317         if(sourceJet && targetJet) {    // duplicate check: slot migth be NULL
318             Double_t targetPt(0);
319             Int_t iSJ(sourceJet->GetNumberOfTracks());
320             Int_t iTJ(targetJet->GetNumberOfTracks());
321             Int_t overlap(0), alreadyFound(0);
322             for(Int_t j(0); j < iSJ; j++) {
323                 alreadyFound = 0;
324                 Int_t idSource((Int_t)sourceJet->TrackAt(j));
325                 for(Int_t k(0); k < iTJ; k++) { // compare all target tracks to the source track
326                     if(idSource == targetJet->TrackAt(k) && alreadyFound == 0) {
327                         overlap++;
328                         alreadyFound++; // avoid possible duplicate matching
329                         AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
330                         if(vp) targetPt += vp->Pt();
331                         continue;
332                     }
333                 }
334             }
335             if((float)overlap/(float)iSJ < fMinFracRecoveredConstituents || targetPt/sourceJet->Pt() < fMinFracRecoveredConstituentPt) {
336                 if(fDebug > 0) printf("  \n > Purging jet, recovered constituents ratio  %i / %i = %.2f <  or pt ratio %.2f / %.2f = %.2f < %.2f", 
337                         overlap, iSJ, (float)overlap/(float)iSJ, targetPt, sourceJet->Pt(), targetPt/sourceJet->Pt(), fMinFracRecoveredConstituentPt);
338                 fMatchedJetContainer[i][0] = 0x0;
339                 fMatchedJetContainer[i][1] = 0x0;
340                 continue;
341             }
342             if(sourceJet->Pt() > 0) {
343                 Double_t sourceRho(0), targetRho(0);
344                 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(sourceJet->Phi(), fSourceRadius)*sourceJet->Area();
345                 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(targetJet->Phi(), fTargetRadius)*targetJet->Area();
346                 fProfFracPtMatched->Fill(sourceJet->Pt()-sourceRho, (targetPt-targetRho) / (sourceJet->Pt()-sourceRho));
347                 fProfFracPtJets->Fill(sourceJet->Pt()-sourceRho, (targetJet->Pt()-targetRho) / (sourceJet->Pt()-sourceRho));
348                 fProfFracNoMatched->Fill(sourceJet->Pt()-sourceRho, (double)overlap / (double)sourceJet->GetNumberOfTracks());
349                 fProfFracNoJets->Fill(sourceJet->Pt()-sourceRho, (double)targetJet->GetNumberOfTracks() / (double)sourceJet->GetNumberOfTracks());
350             }
351             if(fDebug > 0) {
352                 printf("\n > Jet A: %i const\t", iSJ);
353                 printf(" > Jet B %i const\t", iTJ);
354                 printf(" -> OVERLAP: %i tracks <- \n", overlap);
355             }
356         }
357     }
358 }
359 //_____________________________________________________________________________
360 void AliAnalysisTaskJetMatching::GetBijection()
361 {
362     // bijection of source and matched jets, based on closest distance between jets
363     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
364     for(Int_t i(0); i < fNoMatchedJets; i++) {
365         for(Int_t j(i+1); j < fNoMatchedJets; j++) {
366             if(fMatchedJetContainer[i][0] == fMatchedJetContainer[j][0]) {
367                 // found source with two targets, now see which target is closer to the source
368                 if(!(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1] && fMatchedJetContainer[j][0] && fMatchedJetContainer[j][1] )) continue;
369                 Double_t rA(GetR(fMatchedJetContainer[i][0], fMatchedJetContainer[i][1]));      // distance between connection A = i
370                 Double_t rB(GetR(fMatchedJetContainer[j][0], fMatchedJetContainer[j][1]));      // distance between connection B = j
371                 if (rB > rA) {  // jet two is far away, clear it from both target and source list
372                     fMatchedJetContainer[j][0] = 0x0;
373                     fMatchedJetContainer[j][1] = 0x0;
374                 } else {                // jet one is far away, clear it from both target and source list
375                     fMatchedJetContainer[i][0] = fMatchedJetContainer[j][0];    // switch to avoid breaking loop
376                     fMatchedJetContainer[i][1] = fMatchedJetContainer[j][1];    
377                     fMatchedJetContainer[j][0] = 0x0;                           // clear duplicate jet from cache
378                     fMatchedJetContainer[j][1] = 0x0;
379                 }
380                 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
381             }
382         }
383     }
384 }
385 //_____________________________________________________________________________
386 void AliAnalysisTaskJetMatching::PostMatchedJets()
387 {
388     // post matched jets
389     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
390     for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
391         if(fMatchedJetContainer[i][1]) {        // duplicate jets will have NULL value here and are skipped
392             new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
393             p++; 
394         }
395     }
396 }
397 //_____________________________________________________________________________
398 void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
399 {
400     // fill the analysis summary histrogram, saves all relevant analysis settigns
401     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
402     fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
403     fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
404     fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
405     fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
406     fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
407     fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
408     fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
409     fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
410     fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
411     fHistAnalysisSummary->SetBinContent(5, fMatchEta);
412     fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
413     fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
414     fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
415     fHistAnalysisSummary->SetBinContent(7, fMatchR);
416     fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
417     fHistAnalysisSummary->SetBinContent(8, 1.);
418 }
419 //_____________________________________________________________________________
420 void AliAnalysisTaskJetMatching::FillMatchedJetHistograms() 
421 {
422     // fill matched jet histograms
423     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
424     for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
425         AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
426         if(!source) continue;
427         else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
428         Double_t sourceRho(0), targetRho(0);
429         if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
430         fHistSourceJetPt->Fill(source->Pt()-sourceRho);
431         fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
432         for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
433             AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
434             if(target) {
435                 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
436                 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
437                 fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));  
438                 fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
439                 fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
440                 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
441                 if(j==0) {
442                     fHistTargetJetPt->Fill(target->Pt()-targetRho);
443                     fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
444                 }
445             }
446         }
447     }
448     for(Int_t i(0); i < fNoMatchedJets; i++) {
449         if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
450             Double_t sourceRho(0), targetRho(0);
451             if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
452             if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
453             fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
454             fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
455             fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
456             fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
457             fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
458             fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
459             fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
460         }
461     }
462 }
463 //_____________________________________________________________________________
464 void AliAnalysisTaskJetMatching::PrintInfo() const
465 {
466     // print some info 
467     printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
468     printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
469     printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
470     if(fDebug > 3) InputEvent()->GetList()->ls();
471 }
472 //_____________________________________________________________________________
473 void AliAnalysisTaskJetMatching::Terminate(Option_t *)
474 {
475     // terminate
476 }
477 //_____________________________________________________________________________