Add matching task for embedded jets (from Redmer)
authormvl <mvl@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Aug 2013 12:26:12 +0000 (12:26 +0000)
committermvl <mvl@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Aug 2013 12:26:12 +0000 (12:26 +0000)
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskJetMatching.C [new file with mode: 0644]

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
new file mode 100644 (file)
index 0000000..54f3a64
--- /dev/null
@@ -0,0 +1,389 @@
+// 
+// General task to match two sets of jets
+//
+// This task takes two TClonesArray's as input: 
+// [1] fSourceJets - e.g. pythia jets
+// [2] fTargetJets - e.g. a samle containing pythia jets embedded in a pbpb event
+// the task will try to match jets from the source array to the target array, and
+// save the found TARGET jets in a new array 'fMatchedJets'
+//
+// Author: Redmer Alexander Bertens, Utrecht Univeristy, Utrecht, Netherlands
+// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl 
+
+// root includes
+#include <TClonesArray.h>
+#include <TChain.h>
+#include <TMath.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
+// aliroot includes
+#include <AliAnalysisTask.h>
+#include <AliAnalysisManager.h>
+#include <AliLog.h>
+#include <AliVEvent.h>
+#include <AliVParticle.h>
+// emcal jet framework includes
+#include <AliEmcalJet.h>
+#include <AliAnalysisTaskJetMatching.h>
+
+class AliAnalysisTaskJetMatching;
+using namespace std;
+
+ClassImp(AliAnalysisTaskJetMatching)
+
+AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", kTRUE), 
+    fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fUseScaledRho(0), fMatchingScheme(kGeoEtaPhi), fDuplicateJetRecoveryMode(kDoNothing), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(100), fMatchEta(.03), fMatchPhi(.03), fMatchR(.03), fMatchArea(0), fMaxRelEnergyDiff(.1), fMaxAbsEnergyDiff(5) {
+    // default constructor
+    ClearMatchedJetsCache();
+}
+//_____________________________________________________________________________
+AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
+    fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fUseScaledRho(0), fMatchingScheme(kGeoEtaPhi), fDuplicateJetRecoveryMode(kDoNothing), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(100), fMatchEta(.03), fMatchPhi(.03), fMatchR(.03), fMatchArea(0), fMaxRelEnergyDiff(.1), fMaxAbsEnergyDiff(5) {
+    // constructor
+    ClearMatchedJetsCache();
+    DefineInput(0, TChain::Class());
+    DefineOutput(1, TList::Class());
+}
+//_____________________________________________________________________________
+AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
+{
+    // destructor
+    if(fOutputList)             delete fOutputList;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::ExecOnce() 
+{
+    // initialize the anaysis
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(fLocalJetMinEta > -10 && fLocalJetMaxEta > -10) SetJetEtaLimits(fLocalJetMinEta, fLocalJetMaxEta);
+    if(fLocalJetMinPhi > -10 && fLocalJetMaxPhi > -10) SetJetPhiLimits(fLocalJetMinPhi, fLocalJetMaxPhi);
+    // get the stand alone jets from the input event
+    fSourceJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSourceJetsName.Data()));
+    if(!fSourceJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
+    fTargetJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTargetJetsName.Data()));
+    if(!fTargetJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
+    // append the list of matched jets to the event
+    fMatchedJets->Delete();
+    if (!(InputEvent()->FindListObject(fMatchedJetsName))) InputEvent()->AddObject(fMatchedJets);
+    else AliFatal(Form("%s: Object with name %s already in event! Aborting", GetName(), fMatchedJetsName.Data()));
+    FillAnalysisSummaryHistogram();
+    AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
+{
+    // create output objects
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+     // add the matched jets array to the event
+    fMatchedJets = new TClonesArray("AliEmcalJet");
+    fMatchedJets->SetName(fMatchedJetsName);
+    fOutputList = new TList();
+    fOutputList->SetOwner(kTRUE);
+    // add analysis histograms
+    fHistUnsortedCorrelation = BookTH1F("fHistUnsortedCorrelation", "#Delta #varphi unsorted", 50, 0, TMath::Pi());
+    fHistMatchedCorrelation = BookTH1F("fHistMatchedCorrelation", "#Delta #varphi matched", 50, 0, TMath::Pi());
+    fHistSourceJetPt = BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 50, 0, 150);
+    fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 50, 0, 150);
+    fHistMatchedJetPt = BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 50, 0, 150);
+    fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents", 50, 0, 50);
+    fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents", 50, 0, 50);
+    fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents", 50, 0, 50);
+    // the analysis summary histo which stores all the analysis flags is always written to file
+    fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
+    fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
+    fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
+    fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
+    fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
+    fOutputList->Add(fProfQAMatched);
+    fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
+    fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
+    fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
+    fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
+    fOutputList->Add(fProfQA);
+    fOutputList->Sort();
+    PostData(1, fOutputList);
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
+{
+    // book a TH1F and connect it to the output container
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(!fOutputList) return 0x0;
+    TString title(name);
+    title += Form(";%s;[counts]", x);
+    TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
+    histogram->Sumw2();
+    if(append) fOutputList->Add(histogram);
+    return histogram;   
+}
+//_____________________________________________________________________________
+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)
+{
+    // book a TH2F and connect it to the output container
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(!fOutputList) return 0x0;
+    TString title(name);
+    title += Form(";%s;%s", x, y);
+    TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
+    histogram->Sumw2();
+    if(append) fOutputList->Add(histogram);
+    return histogram;   
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskJetMatching::Run()
+{
+    // execute once for each event
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(!(InputEvent() && fSourceJetsName && fTargetJets && IsEventSelected())) return kFALSE;
+    // single event loop
+    switch (fMatchingScheme) {
+        case kGeoEtaPhi : {
+            DoGeometricMatchingEtaPhi();
+        } break;
+        case kGeoR : {
+            DoGeometricMatchingR();
+        } break;
+        case kGeoEtaPhiArea : {
+            DoGeometricMatchingEtaPhi(kTRUE);
+            } break;
+        case kGeoRArea : {
+            DoGeometricMatchingR(kTRUE);
+            } break;
+        case kDeepMatching : {
+            DoGeometricMatchingEtaPhi();
+            } break;
+       default : break;
+    }
+    if(fMatchedJetContainer[1][0]) {       // if matched jets are found, fill some more histograms
+        switch (fDuplicateJetRecoveryMode) {
+            case kDoNothing : break;
+            default : {
+                DuplicateJetRecovery();
+                break;
+            }
+        }
+    }
+    // if required do deep matching, i.e. match constituents in source and target jets 
+    switch (fMatchingScheme) {
+        case kDeepMatching : {
+            DoDeepMatching();
+            break; }
+        default : break;
+    }
+    // stream data to output
+    PostMatchedJets();
+    FillMatchedJetHistograms();
+    if(fDebug > 0) PrintInfo();
+    PostData(1, fOutputList);
+    return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi(Bool_t pairCuts)
+{
+    // do geometric matching
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    fNoMatchedJets = 0; // reset the matched jet counter
+    Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
+    for(Int_t i(0); i < iSource; i++) {
+        AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
+        if(!PassesCuts(sourceJet)) continue;
+        for(Int_t j(0); j < iTarget; j++) {
+            AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
+            if(!PassesCuts(targetJet)) continue;
+            if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta ) && (TMath::Abs(sourceJet->Phi()-targetJet->Phi()) < fMatchPhi)) {
+                if(pairCuts && !PassesCuts(sourceJet, targetJet)) continue;
+                fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
+                fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
+                fNoMatchedJets++;
+                if(fNoMatchedJets > 99) {
+                    AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
+                    return;
+                }
+            }
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::DoGeometricMatchingR(Bool_t pairCuts)
+{
+    // do geometric matching based on shortest path between jet centers
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    fNoMatchedJets = 0; // reset the matched jet counter
+    Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
+    for(Int_t i(0); i < iSource; i++) {
+        AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
+        if(!PassesCuts(sourceJet)) continue;
+        for(Int_t j(0); j < iTarget; j++) {
+            AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
+            if(!PassesCuts(targetJet)) continue;
+            Double_t etaS(sourceJet->Eta()), etaT(targetJet->Eta());
+            Double_t phiS(sourceJet->Phi()), phiT(targetJet->Phi());
+            // if necessary change phase
+            if(TMath::Abs(phiS-phiT) > TMath::Abs(phiS-phiT + TMath::TwoPi())) phiS+=TMath::TwoPi();
+            if(TMath::Abs(phiS-phiT) > TMath::Abs(phiS-phiT - TMath::TwoPi())) phiS-=TMath::TwoPi();
+            if(TMath::Sqrt(TMath::Abs((etaS-etaT)*(etaS-etaT)+(phiS-phiT)*(phiS-phiT)) <= fMatchR)) {
+                if(pairCuts && !PassesCuts(sourceJet, targetJet)) continue;
+                fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
+                fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
+                fNoMatchedJets++;
+                if(fNoMatchedJets > 99) {
+                    AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
+                    return;
+                }
+            }
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::DoDeepMatching()
+{
+    // match constituents, can be VERY slow 
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    if(!fTracks) {
+        AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
+        return; // coverity ...
+    }
+    for(Int_t i(0); i < fNoMatchedJets; i++) {
+        AliEmcalJet* sourceJet = fMatchedJetContainer[0][i];
+        AliEmcalJet* targetJet = fMatchedJetContainer[1][i];
+        if(sourceJet && targetJet) {    // duplicate check: slot migth be NULL
+            Int_t iSJ = sourceJet->GetNumberOfTracks();
+            Int_t iTJ = targetJet->GetNumberOfTracks();
+            for(Int_t j(0); j < iSJ; j++) {     // nested loops over constituets ...
+                AliVParticle* pSJ = sourceJet->TrackAt(j, fTracks);
+                if(!pSJ) continue;              // shouldn't happen
+                for(Int_t k(0); k < iTJ; k++) {
+                    AliVParticle* pTJ = targetJet->TrackAt(k, fTracks);
+                    if(!pTJ) continue;
+                    if(pTJ == pSJ) printf(" > matched track by pointer value \n");
+                    else if(CompareTracks(pSJ, pTJ)) printf(" > hurray, matched source and target constituent ! \n");
+                }
+            }
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::DuplicateJetRecovery()
+{
+    // find target jets that have been matched to a source jet more than once - uses nested loops!
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    Int_t iDuplicateJets(0);            // counter for duplicate jets
+    for(Int_t i(0); i < fNoMatchedJets; i++) {
+        for(Int_t j(i+1); j < fNoMatchedJets; j++) {
+            if(fMatchedJetContainer[i][1] == fMatchedJetContainer[j][1]) {
+                iDuplicateJets++;
+                switch (fDuplicateJetRecoveryMode) {
+                    case kTraceDuplicates : { 
+                        printf(" > found duplicate jet <\n");
+                        break; }
+                    case kRemoveDuplicates : { 
+                         fMatchedJetContainer[j][1] = NULL;
+                         break; }
+                    default : break;
+                }
+            }
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::PostMatchedJets()
+{
+    // post matched jets
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
+        if(fMatchedJetContainer[i][1]) {        // duplicate jets will have NULL value here and are skipped
+            new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
+            p++; 
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
+{
+    // fill the analysis summary histrogram, saves all relevant analysis settigns
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
+    fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
+    fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
+    fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
+    fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fLocalJetMinEta");
+    fHistAnalysisSummary->SetBinContent(5, fLocalJetMinEta);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fLocalJetMaxEta");
+    fHistAnalysisSummary->SetBinContent(6, fLocalJetMaxEta);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fLocalJetMinPhi");
+    fHistAnalysisSummary->SetBinContent(7, fLocalJetMinPhi);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "fLocalJetMaxPhi");
+    fHistAnalysisSummary->SetBinContent(8, fLocalJetMaxPhi);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(9, "fMatchEta");
+    fHistAnalysisSummary->SetBinContent(9, fMatchEta);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(10, "fMatchPhi");
+    fHistAnalysisSummary->SetBinContent(10, fMatchPhi);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(11, "fMatchR");
+    fHistAnalysisSummary->SetBinContent(11, fMatchR);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(12, "fMatchArea");
+    fHistAnalysisSummary->SetBinContent(12, fMatchArea);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(13, "fMaxRelEnergyDiff");
+    fHistAnalysisSummary->SetBinContent(13, fMaxRelEnergyDiff);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(13, "fMaxAbsEnergyDiff");
+    fHistAnalysisSummary->SetBinContent(13, fMaxAbsEnergyDiff);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(14, "iter");
+    fHistAnalysisSummary->SetBinContent(13, 1.);
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::FillMatchedJetHistograms() const
+{
+    // fill matched jet histograms
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
+        AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
+        if(!source) continue;
+        fHistSourceJetPt->Fill(source->Pt());
+        fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
+        for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
+            AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
+            if(target) {
+            fProfQA->Fill(0.5, TMath::Abs(source->Pt()-target->Pt()));
+            fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
+            fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
+                fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
+                if(j==0) {
+                    fHistTargetJetPt->Fill(target->Pt());
+                    fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
+                }
+            }
+        }
+    }
+    for(Int_t i(0); i < fNoMatchedJets; i++) {
+        if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
+            fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
+            fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt());
+            fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt());
+            fProfQAMatched->Fill(0.5, TMath::Abs(fMatchedJetContainer[i][0]->Pt()-fMatchedJetContainer[i][1]->Pt()));
+            fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
+            fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
+        }
+    }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::PrintInfo() const
+{
+    // print some info 
+    printf(" > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
+    printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
+    printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
+    if(fDebug > 3) InputEvent()->GetList()->ls();
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::Terminate(Option_t *)
+{
+    // terminate
+}
+//_____________________________________________________________________________
diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.h
new file mode 100644 (file)
index 0000000..d87b6ad
--- /dev/null
@@ -0,0 +1,143 @@
+#ifndef AliAnalysisTaskJetMatching_H
+#define AliAnalysisTaskJetMatching_H
+
+#include <AliAnalysisTaskEmcalJet.h>
+#include <AliEmcalJet.h>
+#include <AliVTrack.h>
+#include <TClonesArray.h>
+#include <TMath.h>
+#include <TRandom3.h>
+
+class TF1;
+class THF1;
+class THF2;
+class TProfile;
+class AliRhoParameter;
+class AliLocalRhoParameter;
+class TClonesArray;
+
+class AliAnalysisTaskJetMatching : public AliAnalysisTaskEmcalJet
+{
+    public:
+        // enumerators
+        enum matchingSceme      {kGeoEtaPhi, kGeoR, kGeoEtaPhiArea, kGeoRArea, kDeepMatching};
+        enum duplicateRecovery  {kDoNothing, kTraceDuplicates, kRemoveDuplicates}; 
+        enum sourceBKG          {kNoSourceBKG, kSourceRho, kSourceLocalRho};
+        enum targetBKG          {kNoTargetBKG, kTargetRho, kTargetLocalRho};
+        // constructors, destructor
+                                AliAnalysisTaskJetMatching();
+                                AliAnalysisTaskJetMatching(const char *name);
+        virtual                 ~AliAnalysisTaskJetMatching();
+        // setting up the task and technical aspects
+        void                    ExecOnce();
+        virtual void            UserCreateOutputObjects();
+        TH1F*                   BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append = kTRUE);
+        TH2F*                   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 = kTRUE);
+        virtual Bool_t          Run();
+        /* inline */    Double_t PhaseShift(Double_t x) const {  
+            while (x>=TMath::TwoPi())x-=TMath::TwoPi();
+            while (x<0.)x+=TMath::TwoPi();
+            return x; }
+        /* inline */    Double_t PhaseShift(Double_t x, Double_t n) const {
+            x = PhaseShift(x);
+            if(TMath::Nint(n)==2) while (x>TMath::Pi()) x-=TMath::Pi();
+            if(TMath::Nint(n)==3) {
+                if(x>2.*TMath::TwoPi()/n) x = TMath::TwoPi() - x;
+                if(x>TMath::TwoPi()/n) x = TMath::TwoPi()-(x+TMath::TwoPi()/n);
+            }
+            return x; }
+        /* inline */    Bool_t CompareTracks(AliVParticle* a, AliVParticle* b) const {
+//            printf("dEta %.2f \t dPhi %.2f \n", a->Eta()-b->Eta(), a->Phi()-b->Phi());
+            return (TMath::AreEqualAbs(a->Pt(), b->Pt(), 1e-2) && TMath::AreEqualAbs(a->Eta(), b->Eta(), 1e-2) && TMath::AreEqualAbs(a->Phi(), b->Phi(), 1e-2)) ? kTRUE : kFALSE; }
+
+        // setters - setup how to run
+        void                    SetDebugMode(Int_t d)                           {fDebug = d;}
+        void                    SetMatchingScheme(matchingSceme m)              {fMatchingScheme = m;}
+        void                    SetDuplicateRecoveryScheme(duplicateRecovery d) {fDuplicateJetRecoveryMode = d;}
+        void                    SetSourceBKG(sourceBKG b)                       {fSourceBKG = b;}
+        void                    SetTargetBKG(targetBKG b)                       {fTargetBKG = b;}
+        void                    SetSourceJetsName(const char* n)                {fSourceJetsName = n;}
+        void                    SetTargetJetsName(const char* n)                {fTargetJetsName = n; }
+        void                    SetMatchedJetsName(const char* n)               {fMatchedJetsName = n;}
+        void                    SetMatchEta(Float_t f)                          {fMatchEta = f;}
+        void                    SetMatchPhi(Float_t f)                          {fMatchPhi = f;}
+        void                    SetMatchR(Float_t f)                            {fMatchR = f;}
+        void                    SetMatchRelArea(Float_t f)                      {fMatchArea = f;}
+        void                    SetMaxRelEnergyDiff(Float_t f)                  {fMaxRelEnergyDiff = f;}
+        void                    SetMaxAbsEnergyDiff(Float_t f)                  {fMaxAbsEnergyDiff = f;}
+        // methods
+        void                    DoGeometricMatchingEtaPhi(Bool_t pairCuts = kFALSE);
+        void                    DoGeometricMatchingR(Bool_t pairCuts = kFALSE);
+        void                    DoDeepMatching();
+        void                    DuplicateJetRecovery();
+        void                    PostMatchedJets();
+        // fill histogrmas
+        void                    FillAnalysisSummaryHistogram() const;
+        void                    FillMatchedJetHistograms() const;
+        // setters - analysis details
+        /* inline */    Bool_t PassesCuts(const AliVTrack* track) const {
+            if(!track) return kFALSE;
+            return (track->Pt() < fTrackPtCut || track->Eta() < fTrackMinEta || track->Eta() > fTrackMaxEta || track->Phi() < fTrackMinPhi || track->Phi() > fTrackMaxPhi) ? kFALSE : kTRUE; }
+        /* inline */    Bool_t PassesCuts(AliEmcalJet* jet) const {
+            return (jet) ? kTRUE : kFALSE; }
+        /* inline */    Bool_t PassesCuts(AliEmcalJet* a, AliEmcalJet* b) const {
+            if (fMatchArea > 0) { return (TMath::Abs(a->Area()/b->Area()) < fMatchArea) ? kTRUE : kFALSE; }
+            if (fMaxRelEnergyDiff > 0) { return (TMath::Abs(a->E()/b->E()) > fMaxRelEnergyDiff) ? kTRUE : kFALSE; }
+            if (fMaxAbsEnergyDiff > 0) { return (TMath::Abs(a->E()-b->E()) > fMaxAbsEnergyDiff) ? kTRUE : kFALSE; }
+            return kTRUE; }
+        /* inline */    void ClearMatchedJetsCache() {
+            for (Int_t i(0); i < fNoMatchedJets; i++) {
+                fMatchedJetContainer[i][0] = 0x0; fMatchedJetContainer[i][1] = 0x0; }
+            fNoMatchedJets = 0;
+        }
+        void                    PrintInfo() const;
+        virtual void            Terminate(Option_t* option);
+
+    private: 
+        Int_t                   fDebug;                 // debug level (0 none, 1 fcn calls, 2 verbose)
+        TClonesArray*           fSourceJets;            //! array with source jets
+        TString                 fSourceJetsName;        // name of array with source jets
+        TClonesArray*           fTargetJets;            //! array with target jets
+        TString                 fTargetJetsName;        // name of array with target jets
+        TClonesArray*           fMatchedJets;           //! final list of matched jets which is added to event
+        TString                 fMatchedJetsName;       // name of list of matched jets
+        Bool_t                  fUseScaledRho;          // use scaled rho
+        matchingSceme           fMatchingScheme;        // select your favorite matching algorithm
+        duplicateRecovery       fDuplicateJetRecoveryMode;      // what to do with duplicate matches
+        sourceBKG               fSourceBKG;             // subtracted background of source jets
+        targetBKG               fTargetBKG;             // subtracted background of target jets
+        // additional jet cuts (most are inherited)
+        Float_t                 fLocalJetMinEta;        // local eta cut for jets
+        Float_t                 fLocalJetMaxEta;        // local eta cut for jets
+        Float_t                 fLocalJetMinPhi;        // local phi cut for jets
+        Float_t                 fLocalJetMaxPhi;        // local phi cut for jets
+        // transient object pointers
+        TList*                  fOutputList;            //! output list
+        TH1F*                   fHistUnsortedCorrelation;       //! dphi correlation of unsorted jets
+        TH1F*                   fHistMatchedCorrelation;        //! dphi correlation of matched jets
+        TH1F*                   fHistSourceJetPt;       //! pt of source jets
+        TH1F*                   fHistTargetJetPt;       //! pt of target jets
+        TH1F*                   fHistMatchedJetPt;      //! pt of matched jets
+        TH1F*                   fHistNoConstSourceJet;  //! number of constituents of source jet
+        TH1F*                   fHistNoConstTargetJet;  //! number of constituents of target jet
+        TH1F*                   fHistNoConstMatchJet;   //! number of constituents of matched jets
+        TH1F*                   fHistAnalysisSummary;   //! flags
+        TProfile*               fProfQAMatched;         //! QA spreads of matched jets
+        TProfile*               fProfQA;                //! QA spreads of source and target jets
+        Int_t                   fNoMatchedJets;         //! number of matched jets
+        AliEmcalJet*            fMatchedJetContainer[100][2];   //! pointers to matched jets
+        // geometric matching parameters
+        Float_t                 fMatchEta;              // max eta distance between centers of matched jets
+        Float_t                 fMatchPhi;              // max phi distance between centers of matched jets
+        Float_t                 fMatchR;                // max distance between matched jets
+        Float_t                 fMatchArea;             // max relative area mismatch between matched jets
+        Float_t                 fMaxRelEnergyDiff;      // max relative energy difference between matched jets
+        Float_t                 fMaxAbsEnergyDiff;      // max absolute energy difference between matched jets
+
+        AliAnalysisTaskJetMatching(const AliAnalysisTaskJetMatching&);                  // not implemented
+        AliAnalysisTaskJetMatching& operator=(const AliAnalysisTaskJetMatching&);       // not implemented
+
+        ClassDef(AliAnalysisTaskJetMatching, 1);
+};
+
+#endif
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskJetMatching.C b/PWGJE/EMCALJetTasks/macros/AddTaskJetMatching.C
new file mode 100644 (file)
index 0000000..dbce90a
--- /dev/null
@@ -0,0 +1,36 @@
+/* AddTask macro for class AddTaskJetMatching.C
+ * Redmer Alexander Bertens, rbertens@cern.ch
+ * Utrecht University, Utrecht, Netherlands             */
+
+AliAnalysisTaskJetMatching* AddTaskJetMatching(
+        const char* sourceJets  = "SourceJets", // source jets
+        const char* targetJets  = "TargetJets", // target jets
+        const char* matchedJets = "MatchedJets",// matched jets
+        UInt_t matchingScheme   = AliAnalysisTaskJetMatching::kDeepMatching,
+        UInt_t duplicateRecovery= AliAnalysisTaskJetMatching::kTraceDuplicates,
+        const char *name        = "AliAnalysisTaskJetMatching",
+  )
+{ 
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)                             return 0x0;
+  if (!mgr->GetInputEventHandler())     return 0x0;
+
+  AliAnalysisTaskJetMatching* jetTask = new AliAnalysisTaskJetMatching(name);
+  jetTask->SetDebugMode(-1);
+  jetTask->SetSourceJetsName(sourceJets);
+  jetTask->SetTargetJetsName(targetJets);
+  jetTask->SetMatchedJetsName(matchedJets);
+  jetTask->SetMatchingScheme(matchingScheme);
+  jetTask->SetDuplicateRecoveryScheme(duplicateRecovery);
+  mgr->AddTask(jetTask);
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer()  ;
+  TString contname(name);
+  contname += "_histos";
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(), 
+                                                           TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                           Form("%s", AliAnalysisManager::GetCommonFileName()));
+  mgr->ConnectInput  (jetTask, 0,  cinput1 );
+  mgr->ConnectOutput (jetTask, 1, coutput1 );
+  return jetTask;
+}