PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx

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