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