]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
bug-fix: rotation of sub-leading jet in di-jet
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetMatching.cxx
index 0851aa0f38c525b06c485d4fad323be695999883..0e3daebb0b949eafeb6ef0719e61f17ce8dca20b 100644 (file)
@@ -1,3 +1,18 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
 // 
 // General task to match two sets of jets
 //
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TProfile.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TKey.h>
+#include <TSystem.h>
 // aliroot includes
 #include <AliAnalysisTask.h>
 #include <AliAnalysisManager.h>
@@ -60,14 +79,14 @@ using namespace std;
 
 ClassImp(AliAnalysisTaskJetMatching)
 
-AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJetDev("AliAnalysisTaskJetMatching", kTRUE), 
-    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) {
+AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", kTRUE), 
+    fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fSourceRho(0), fSourceRhoName(0), fTargetRho(0), fTargetRhoName(0), fUseScaledRho(0), fSourceRadius(0.3), fTargetRadius(0.3), fMatchingScheme(kGeoEtaPhi), fUseEmcalBaseJetCuts(kFALSE), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistMatchedSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistSourceMatchedJetPt(0), fHistDetectorResponseProb(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fProfFracPtMatched(0), fProfFracPtJets(0), fProfFracNoMatched(0), fProfFracNoJets(0), fHistDiJet(0), fHistDiJetLeadingJet(0), fHistDiJetDPhi(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(200), fMatchEta(.3), fMatchPhi(.3), fMatchR(.08), fDoDetectorResponse(kFALSE), fMatchConstituents(kTRUE), fMinFracRecoveredConstituents(.5), fMinFracRecoveredConstituentPt(0.5), fGetBijection(kTRUE), fh1Trials(0x0), fh1AvgTrials(0x0), fh1Xsec(0x0), fAvgTrials(0) {
     // default constructor
     ClearMatchedJetsCache();
 }
 //_____________________________________________________________________________
-AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJetDev(name, kTRUE),
-    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) {
+AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
+    fDebug(0), fSourceJets(0), fSourceJetsName(0), fTargetJets(0), fTargetJetsName(0), fMatchedJets(0), fMatchedJetsName(GetName()), fSourceRho(0), fSourceRhoName(0), fTargetRho(0), fTargetRhoName(0), fUseScaledRho(0), fSourceRadius(0.3), fTargetRadius(0.3), fMatchingScheme(kGeoEtaPhi), fUseEmcalBaseJetCuts(kFALSE), fSourceBKG(kNoSourceBKG), fTargetBKG(kNoTargetBKG), fOutputList(0), fHistUnsortedCorrelation(0), fHistMatchedCorrelation(0), fHistSourceJetPt(0), fHistMatchedSourceJetPt(0), fHistTargetJetPt(0), fHistMatchedJetPt(0), fHistSourceMatchedJetPt(0), fHistDetectorResponseProb(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fProfFracPtMatched(0), fProfFracPtJets(0), fProfFracNoMatched(0), fProfFracNoJets(0), fHistDiJet(0), fHistDiJetLeadingJet(0), fHistDiJetDPhi(0), fHistAnalysisSummary(0), fProfQAMatched(0), fProfQA(0), fNoMatchedJets(200), fMatchEta(.3), fMatchPhi(.3), fMatchR(.08), fDoDetectorResponse(kFALSE), fMatchConstituents(kTRUE), fMinFracRecoveredConstituents(0.5), fMinFracRecoveredConstituentPt(0.5), fGetBijection(kTRUE), fh1Trials(0x0), fh1AvgTrials(0x0), fh1Xsec(0x0), fAvgTrials(0) {
     // constructor
     ClearMatchedJetsCache();
     DefineInput(0, TChain::Class());
@@ -108,7 +127,7 @@ void AliAnalysisTaskJetMatching::ExecOnce()
         } break;
         default : break;
     }
-    AliAnalysisTaskEmcalJetDev::ExecOnce(); // init base class
+    AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
     if(fDoDetectorResponse) fMatchConstituents = kFALSE;
 }
 //_____________________________________________________________________________
@@ -128,14 +147,26 @@ void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
     fHistMatchedSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistMatchedParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistMatchedSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
     fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
     fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
-    fHistSourceMatchedJetPt = (fDoDetectorResponse) ? BookTH2F("fHistPartDetJetPt", "particle level jet p_{t}^{gen} [GeV/c]", "detector level jet p_{t}^{rec} [GeV/c]", 150, -50, 250, 150, -50, 250) : BookTH2F("fHistSourceMatchedJetPt", "source jet p_{t} [GeV/c]", "matched jet p_{t} [GeV/c]", 150, -50, 250, 150, -50, 250);
-    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);
+    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);
+    if(fDoDetectorResponse) {
+        fHistDetectorResponseProb = BookTH2F("fHistDetectorResponseProb", "(p_{t}^{det} - p_{t}^{part})/p_{t}^{part}", "p_{t}^{part}", 100, -1.5, 1., 20, 0, 200);
+        fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+        fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+        fOutputList->Add(fh1Xsec);
+        fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
+        fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+        fOutputList->Add(fh1Trials);
+        fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
+        fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
+        fOutputList->Add(fh1AvgTrials);
+    }
+    fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents source jets", 50, 0, 50);
+    fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents target jets", 50, 0, 50);
+    fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents matched jets", 50, 0, 50);
     if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
-        fProfFracPtMatched = new TProfile("fProfFracPtMatched", "fProfFracPtMatched", 15, -50, 250);
+        fProfFracPtMatched = new TProfile("fProfFracPtMatched", "recovered target p_{T} / source p_{T}", 15, -50, 250);
         fOutputList->Add(fProfFracPtMatched);
-        fProfFracNoMatched = new TProfile("fProfFracNoMatched", "fProfFracNoMatched", 15, -50, 250);
+        fProfFracNoMatched = new TProfile("fProfFracNoMatched", "recovered target constituents / source constituents", 15, -50, 250);
         fOutputList->Add(fProfFracNoMatched);
     }
     // the analysis summary histo which stores all the analysis flags is always written to file
@@ -150,10 +181,18 @@ void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
     fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
     fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
     fOutputList->Add(fProfQA);
-    fProfFracPtJets = new TProfile("fProfFracPtJets", "fProfFracPtJets", 15, -50, 250);
+    fProfFracPtJets = new TProfile("fProfFracPtJets", "source p_{T} / target p_{T}", 15, -50, 250);
     fOutputList->Add(fProfFracPtJets);
-    fProfFracNoJets = new TProfile("fProfFracNoJets", "fProfFracNoJets", 15, -50, 250);
+    fProfFracNoJets = new TProfile("fProfFracNoJets", "source constituents / target constituents", 15, -50, 250);
     fOutputList->Add(fProfFracNoJets);
+    switch (fMatchingScheme) {
+        case kDiJet : {
+            fHistDiJet = BookTH2F("fHistDiJet", "matched di-jet #varphi", "matched di-jet #eta", 100, 0., TMath::TwoPi(), 100, -5., 5.);
+            fHistDiJetLeadingJet = BookTH2F("fHistDiJetLeadingJet", "leading jet #varphi", "leadingd jet #eta", 100, 0., TMath::TwoPi(), 100, -5., 5.);
+            fHistDiJetDPhi = BookTH1F("fHistDiJetDPhi", "leading jet #varphi - (matched jet #varphi - #pi)", 100, -1.*TMath::Pi(), TMath::Pi());
+        } break;
+        default : break;
+    }
     fOutputList->Sort();
     PostData(1, fOutputList);
 }
@@ -184,19 +223,77 @@ TH2F* AliAnalysisTaskJetMatching::BookTH2F(const char* name, const char* x, cons
     return histogram;   
 }
 //_____________________________________________________________________________
+Bool_t AliAnalysisTaskJetMatching::Notify()
+{
+    // for each file get the number of trials and pythia cross section
+    // see  $ALICE_ROOT/PWGJE/AliAnalysisTaskJetProperties.cxx
+    //      $ALICE_ROOT/PWG/Tools/AliAnalysisHelperJetTasks.cxx
+    // this function is implenented here temporarily to avoid introducing a dependency
+    // later on this could just be a call to a static helper function
+    if(!fDoDetectorResponse) return kTRUE;
+    Float_t xsection(0), ftrials(1);
+    fAvgTrials = ftrials;
+    TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+    if(tree) {
+        TFile *curfile = tree->GetCurrentFile();
+        if (!curfile) return kFALSE;
+        TString file(curfile->GetName()); 
+        if(file.Contains("root_archive.zip#")) file.Replace(file.Index("#",1,file.Index("root_archive",12,0,TString::kExact),TString::kExact)+1,file.Index(".root",5,TString::kExact)-file.Index("root_archive",12,0,TString::kExact),"");
+        else file.ReplaceAll(gSystem->BaseName(file.Data()),"");
+        TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
+        if(!fxsec) {
+            fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
+            if(fxsec) {
+                TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
+                if(key) {
+                      TList *list = dynamic_cast<TList*>(key->ReadObj());
+                      if(list) {
+                          xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+                          ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+                      }
+                }
+                fxsec->Close();
+            }
+        } else {
+            TTree *xtree = (TTree*)fxsec->Get("Xsection");
+            if(xtree) {
+                UInt_t _ftrials  = 0;
+                Double_t _xsection  = 0;
+                xtree->SetBranchAddress("xsection",&_xsection);
+                xtree->SetBranchAddress("ntrials",&_ftrials);
+                xtree->GetEntry(0);
+                ftrials = _ftrials;
+                xsection = _xsection;
+            }
+            fxsec->Close();
+        }
+        fh1Xsec->Fill("<#sigma>",xsection);
+        Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+        if(ftrials >= nEntries && nEntries > 0.) fAvgTrials = ftrials/nEntries;
+        fh1Trials->Fill("#sum{ntrials}",ftrials); 
+    }  
+    return kTRUE;
+}
+//_____________________________________________________________________________
 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() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
+    if(fh1AvgTrials) fh1AvgTrials->Fill("#sum{avg ntrials}", fAvgTrials);
     // step one: do geometric matching 
     switch (fMatchingScheme) {
+        // FIXME having separate dedicated functions is not necessary and historical
+        // cluttered code - should be cleaned up and merged into one
         case kGeoEtaPhi : {
             DoGeometricMatchingEtaPhi();
         } break;
         case kGeoR : {
             DoGeometricMatchingR();
         } break;
+        case kDiJet : {
+            DoDiJetMatching();
+        } break;
         default : break;
     }
     // break if no jet was matched
@@ -221,11 +318,10 @@ void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
     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;
-        if(fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
+        if(!PassesCuts(sourceJet, 0)) continue;
         for(Int_t j(0); j < iTarget; j++) {
             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
-            if(!PassesCuts(targetJet)) continue;
+            if(!PassesCuts(targetJet, 1)) continue;
             if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
             if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
                 Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
@@ -234,11 +330,10 @@ void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
                 if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {       // accept the jets as matching 
                 Bool_t isBestMatch(kTRUE);
                     if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
-                        if(fDebug > 0) printf(" > Entering first bbijection test \n");
+                        if(fDebug > 0) printf(" > Entering first bijection test \n");
                         for(Int_t k(i); k < iSource; k++) {
                             AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
-                            if(!PassesCuts(candidateSourceJet)) continue;
-                            if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
+                            if(PassesCuts(candidateSourceJet, 0)) continue;
                             if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
                             if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
                                 isBestMatch = kFALSE;
@@ -270,11 +365,10 @@ void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
     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;
-        else if (fUseEmcalBaseJetCuts && !AcceptJet(sourceJet, 0)) continue;
+        if(!PassesCuts(sourceJet, 0)) continue;
         for(Int_t j(0); j < iTarget; j++) {
             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
-            if(!PassesCuts(targetJet)) continue;
+            if(!PassesCuts(targetJet, 1)) continue;
             else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
             if(GetR(sourceJet, targetJet) <= fMatchR) {
                 Bool_t isBestMatch(kTRUE);
@@ -282,8 +376,7 @@ void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
                     if(fDebug > 0) printf(" > Entering first bijection test \n");
                     for(Int_t k(i); k < iSource; k++) {
                         AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
-                        if(!PassesCuts(candidateSourceJet)) continue;
-                        if(fUseEmcalBaseJetCuts && !AcceptJet(candidateSourceJet, 0)) continue;
+                        if(!PassesCuts(candidateSourceJet, 0)) continue;
                         if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
                         if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
                             isBestMatch = kFALSE;
@@ -306,6 +399,38 @@ void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
     }
 }
 //_____________________________________________________________________________
+void AliAnalysisTaskJetMatching::DoDiJetMatching()
+{
+    // match dijets. this is in a sense a 'special' mode of the task as both jet supplied jet arrays 
+    // (target and source) are the same jet collection, matching will be performed on distribution in
+    // azimuth
+    // no ouptut array is produced, dedicated dijet histo's are filled in this function instead
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    // get the leading jet in a given acceptance (TPC is default)
+    Int_t leadingJetIndex(-1), subLeadingJetIndex(-1);
+    // retrieve the leading jet, leadingJetIndex points to the leading jet
+    AliEmcalJet* leadingJet(GetLeadingJet(fSourceJets, leadingJetIndex));
+    if(!leadingJet) return;
+    // fill phi and eta of leading jet (should always be in selected acceptance)
+    fHistDiJetLeadingJet->Fill(leadingJet->Phi(), leadingJet->Eta());
+    Double_t sourcePhi(leadingJet->Phi()), targetPhi(-1);
+    // get the sub-leading jet - faster when leading jet is also provided
+    AliEmcalJet* subLeadingJet(GetSubLeadingJet(fSourceJets, leadingJetIndex, subLeadingJetIndex));
+    if(!subLeadingJet) return;
+    else { // check if the sub-leading jet is actually the away-side jet
+        targetPhi = subLeadingJet->Phi() + TMath::Pi();
+        // rotate jets to common phase
+        if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
+        if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
+        if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi+TMath::TwoPi())) targetPhi+=TMath::TwoPi();
+        if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi-TMath::TwoPi())) targetPhi-=TMath::TwoPi();
+        if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {
+            fHistDiJet->Fill(subLeadingJet->Phi(), subLeadingJet->Eta());
+            fHistDiJetDPhi->Fill(sourcePhi-targetPhi);
+        }
+    }
+}
+//_____________________________________________________________________________
 void AliAnalysisTaskJetMatching::DoConstituentMatching()
 {
     // match constituents within matched jets 
@@ -461,15 +586,58 @@ void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
             fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
             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()));
+            
             fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
             if(fDoDetectorResponse) {
                 fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
                 fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
+                fHistDetectorResponseProb->Fill((fMatchedJetContainer[i][1]->Pt()-fMatchedJetContainer[i][0]->Pt())/fMatchedJetContainer[i][0]->Pt(), fMatchedJetContainer[i][0]->Pt());
             }
         }
     }
 }
 //_____________________________________________________________________________
+AliEmcalJet* AliAnalysisTaskJetMatching::GetLeadingJet(TClonesArray* source, Int_t &leadingJetIndex, Double_t etaMin, Double_t etaMax) 
+{
+    // return the leading jet within given acceptance
+    Int_t iJets(source->GetEntriesFast());
+    Double_t pt(0);
+    AliEmcalJet* leadingJet(0x0);
+    for(Int_t i(0); i < iJets; i++) {
+        AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
+        if(jet->Eta() < etaMin || jet->Eta() > etaMax) continue;
+        if(jet->Pt() > pt) {
+           leadingJet = jet;
+           pt = leadingJet->Pt();
+           leadingJetIndex = i;
+        }
+    }
+    return leadingJet;
+}
+//_____________________________________________________________________________
+AliEmcalJet* AliAnalysisTaskJetMatching::GetSubLeadingJet(TClonesArray* source, Int_t leadingJetIndex, Int_t &subLeadingJetIndex, Double_t etaMin, Double_t etaMax) 
+{
+    // return the leading jet within given acceptance
+    // same as GetLeadingJet() but skips the leading jet (so returned jet is 
+    // sub-leading by design)
+    Int_t iJets(source->GetEntriesFast());
+    // if the leading jet isn't given, retrieve it
+    if(leadingJetIndex < 0) GetLeadingJet(source, leadingJetIndex, etaMin, etaMax);
+    Double_t pt(0);
+    AliEmcalJet* leadingJet(0x0);
+    for(Int_t i(0); i < iJets; i++) {
+        if(i == leadingJetIndex) continue;
+        AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
+        if(jet->Eta() < etaMin || jet->Eta() > etaMax) continue;
+        if(jet->Pt() > pt) {
+           leadingJet = jet;
+           pt = leadingJet->Pt();
+           subLeadingJetIndex = i;
+        }
+    }
+    return leadingJet;
+}
+//_____________________________________________________________________________
 void AliAnalysisTaskJetMatching::PrintInfo() const
 {
     // print some info