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 <TH3.h>
  60 #include <TH1F.h>
  61 #include <TH2F.h>
  62 #include <TH3F.h>
  63 #include <TProfile.h>
  64 #include <TFile.h>
  65 #include <TTree.h>
  66 #include <TKey.h>
  67 #include <TSystem.h>
  68 // aliroot includes
  69 #include <AliAnalysisTask.h>
  70 #include <AliAnalysisManager.h>
  71 #include <AliLog.h>
  72 #include <AliVEvent.h>
  73 #include <AliVParticle.h>
  74 // emcal jet framework includes
  75 #include <AliEmcalJet.h>
  76 #include <AliAnalysisTaskJetMatching.h>
  77 #include <AliLocalRhoParameter.h>
  78
  79 class AliAnalysisTaskJetMatching;
  80 using namespace std;
  81
  82 ClassImp(AliAnalysisTaskJetMatching)
  83
  84 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching() : AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMatching", 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), fHistDetectorResponseProb(0), fHistNoConstSourceJet(0), fHistNoConstTargetJet(0), fHistNoConstMatchJet(0), fProfFracPtMatched(0), fProfFracPtJets(0), fProfFracNoMatched(0), fProfFracNoJets(0), fHistDiJet(0), fHistDiJetLeadingJet(0), fHistDiJetDPhi(0), fHistDiJetDPt(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) {
  86     // default constructor
  87     ClearMatchedJetsCache();
  88 }
  89 //_____________________________________________________________________________
  90 AliAnalysisTaskJetMatching::AliAnalysisTaskJetMatching(const char* name) : AliAnalysisTaskEmcalJet(name, kTRUE),
  91     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), fHistDiJetDPt(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) {
  92     // constructor
  93     ClearMatchedJetsCache();
  94     DefineInput(0, TChain::Class());
  95     DefineOutput(1, TList::Class());
  96 }
  97 //_____________________________________________________________________________
  98 AliAnalysisTaskJetMatching::~AliAnalysisTaskJetMatching()
  99 {
 100     // destructor
 101     if(fOutputList)             delete fOutputList;
 102 }
 103 //_____________________________________________________________________________
 104 void AliAnalysisTaskJetMatching::ExecOnce()
 105 {
 106     // initialize the anaysis
 107     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 108     // get the stand alone jets from the input event
 109     fSourceJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSourceJetsName.Data()));
 110     if(!fSourceJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
 111     fTargetJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTargetJetsName.Data()));
 112     if(!fTargetJets) AliFatal(Form("%s: Container with name %s not found. Aborting", GetName(), fSourceJetsName.Data()));
 113     // append the list of matched jets to the event
 114     fMatchedJets->Delete();
 115     if (!(InputEvent()->FindListObject(fMatchedJetsName))) InputEvent()->AddObject(fMatchedJets);
 116     else AliFatal(Form("%s: Object with name %s already in event! Aborting", GetName(), fMatchedJetsName.Data()));
 117     FillAnalysisSummaryHistogram();
 118     switch (fSourceBKG) {
 119         case kSourceLocalRho : {
 120             fSourceRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fSourceRhoName));
 121             if(!fSourceRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fSourceRhoName.Data()));
 122         } break;
 123         default : break;
 124     }
 125     switch (fTargetBKG) {
 126         case kTargetLocalRho : {
 127             fTargetRho = dynamic_cast<AliLocalRhoParameter*>(InputEvent()->FindListObject(fTargetRhoName));
 128             if(!fTargetRho) AliFatal(Form("%s: Object with name %s requested but not found! Aborting", GetName(), fTargetRhoName.Data()));
 129         } break;
 130         default : break;
 131     }
 132     AliAnalysisTaskEmcalJet::ExecOnce(); // init base class
 133     if(fDoDetectorResponse) fMatchConstituents = kFALSE;
 134 }
 135 //_____________________________________________________________________________
 136 void AliAnalysisTaskJetMatching::UserCreateOutputObjects()
 137 {
 138     // create output objects
 139     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 140      // add the matched jets array to the event
 141     fMatchedJets = new TClonesArray("AliEmcalJet");
 142     fMatchedJets->SetName(fMatchedJetsName);
 143     fOutputList = new TList();
 144     fOutputList->SetOwner(kTRUE);
 145     // add analysis histograms
 146     fHistUnsortedCorrelation = BookTH1F("fHistUnsortedCorrelation", "#Delta #varphi unsorted", 50, 0, TMath::Pi());
 147     fHistMatchedCorrelation = BookTH1F("fHistMatchedCorrelation", "#Delta #varphi matched", 50, 0, TMath::Pi());
 148     fHistSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) :  BookTH1F("fHistSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
 149     fHistMatchedSourceJetPt = (fDoDetectorResponse) ? BookTH1F("fHistMatchedParticleLevelJetPt", "p_{t}^{gen} [GeV/c]", 150, -50, 250) : BookTH1F("fHistMatchedSourceJetPt", "p_{t} [GeV/c]", 150, -50, 250);
 150     fHistTargetJetPt = BookTH1F("fHistTargetJetPt", "p_{t} [GeV/c]", 150, -50, 250);
 151     fHistMatchedJetPt = (fDoDetectorResponse) ? BookTH1F("fHistDetectorLevelJet", "p_{t}^{rec}", 150, -50, 250) : BookTH1F("fHistMatchedJetPt", "p_{t} [GeV/c]", 150, -50, 250);
 152     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);
 153     if(fDoDetectorResponse) {
 154         fHistDetectorResponseProb = BookTH2F("fHistDetectorResponseProb", "(p_{t}^{det} - p_{t}^{part})/p_{t}^{part}", "p_{t}^{part}", 100, -1.5, 1., 20, 0, 200);
 155         fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 156         fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 157         fOutputList->Add(fh1Xsec);
 158         fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 159         fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 160         fOutputList->Add(fh1Trials);
 161         fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
 162         fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 163         fOutputList->Add(fh1AvgTrials);
 164     }
 165     fHistNoConstSourceJet = BookTH1F("fHistNoConstSourceJet", "number of constituents source jets", 50, 0, 50);
 166     fHistNoConstTargetJet = BookTH1F("fHistNoConstTargetJet", "number of constituents target jets", 50, 0, 50);
 167     fHistNoConstMatchJet = BookTH1F("fHistNoConstMatchJet", "number of constituents matched jets", 50, 0, 50);
 168     if(!fDoDetectorResponse) { // these observables cannot be measured in current detector response mode
 169         fProfFracPtMatched = new TProfile("fProfFracPtMatched", "recovered target p_{T} / source p_{T}", 15, -50, 250);
 170         fOutputList->Add(fProfFracPtMatched);
 171         fProfFracNoMatched = new TProfile("fProfFracNoMatched", "recovered target constituents / source constituents", 15, -50, 250);
 172         fOutputList->Add(fProfFracNoMatched);
 173     }
 174     // the analysis summary histo which stores all the analysis flags is always written to file
 175     fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 51, -0.5, 15.5);
 176     fProfQAMatched = new TProfile("fProfQAMatched", "fProfQAMatched", 3, 0, 3);
 177     fProfQAMatched->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
 178     fProfQAMatched->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
 179     fProfQAMatched->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
 180     fOutputList->Add(fProfQAMatched);
 181     fProfQA = new TProfile("fProfQA", "fProfQA", 3, 0, 3);
 182     fProfQA->GetXaxis()->SetBinLabel(1, "<#delta p_{t} >");
 183     fProfQA->GetXaxis()->SetBinLabel(2, "<#delta #eta>");
 184     fProfQA->GetXaxis()->SetBinLabel(3, "<#delta #varphi>");
 185     fOutputList->Add(fProfQA);
 186     fProfFracPtJets = new TProfile("fProfFracPtJets", "source p_{T} / target p_{T}", 15, -50, 250);
 187     fOutputList->Add(fProfFracPtJets);
 188     fProfFracNoJets = new TProfile("fProfFracNoJets", "source constituents / target constituents", 15, -50, 250);
 189     fOutputList->Add(fProfFracNoJets);
 190     switch (fMatchingScheme) {
 191         case kDiJet : {
 192             fHistDiJet = BookTH3F("fHistDiJet", "matched di-jet #varphi", "matched di-jet #eta", "leading jet p_{T} (GeV/c)", 100, 0., TMath::TwoPi(), 100, -5., 5., 100, 0, 200);
 193             fHistDiJetLeadingJet = BookTH3F("fHistDiJetLeadingJet", "leading jet #varphi", "leadingd jet #eta", "leading jet p_{T} (GeV/c)", 100, 0., TMath::TwoPi(), 100, -5., 5., 100, 0, 200);
 194             fHistDiJetDPhi = BookTH2F("fHistDiJetDPhi", "leading jet #varphi - (matched jet #varphi - #pi)", "leading jet p_{T} (GeV/c)", 100, -1.*TMath::Pi(), TMath::Pi(), 100, 0, 200);
 195             fHistDiJetDPt = BookTH2F("fHistDiJetDPt", "leading jet p_{T} - sub leading jet p_{T} (GeV/c)", "leading jet p_{T} (GeV/c)", 100, -25, 25, 100, 0, 200);
 196         } break;
 197         default : break;
 198     }
 199     fOutputList->Sort();
 200     PostData(1, fOutputList);
 201 }
 202 //_____________________________________________________________________________
 203 TH1F* AliAnalysisTaskJetMatching::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Bool_t append)
 204 {
 205     // book a TH1F and connect it to the output container
 206     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 207     if(append && !fOutputList) return 0x0;
 208     TString title(name);
 209     title += Form(";%s;[counts]", x);
 210     TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
 211     histogram->Sumw2();
 212     if(append) fOutputList->Add(histogram);
 213     return histogram;
 214 }
 215 //_____________________________________________________________________________
 216 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)
 217 {
 218     // book a TH2F and connect it to the output container
 219     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 220     if(append && !fOutputList) return 0x0;
 221     TString title(name);
 222     title += Form(";%s;%s", x, y);
 223     TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
 224     histogram->Sumw2();
 225     if(append) fOutputList->Add(histogram);
 226     return histogram;
 227 }
 228 //_____________________________________________________________________________
 229 TH3F* AliAnalysisTaskJetMatching::BookTH3F(const char* name, const char* x, const char* y, const char* z, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t binsz, Double_t minz, Double_t maxz, Bool_t append)
 230 {
 231     // book a TH2F and connect it to the output container
 232     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 233     if(append && !fOutputList) return 0x0;
 234     TString title(name);
 235     title += Form(";%s;%s;%s", x, y, z);
 236     TH3F* histogram = new TH3F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy, binsz, minz, maxz);
 237     histogram->Sumw2();
 238     if(append) fOutputList->Add(histogram);
 239     return histogram;
 240 }
 241 //_____________________________________________________________________________
 242 Bool_t AliAnalysisTaskJetMatching::Notify()
 243 {
 244     // for each file get the number of trials and pythia cross section
 245     // see  $ALICE_ROOT/PWGJE/AliAnalysisTaskJetProperties.cxx
 246     //      $ALICE_ROOT/PWG/Tools/AliAnalysisHelperJetTasks.cxx
 247     // this function is implenented here temporarily to avoid introducing a dependency
 248     // later on this could just be a call to a static helper function
 249     if(!fDoDetectorResponse) return kTRUE;
 250     Float_t xsection(0), ftrials(1);
 251     fAvgTrials = ftrials;
 252     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 253     if(tree) {
 254         TFile *curfile = tree->GetCurrentFile();
 255         if (!curfile) return kFALSE;
 256         TString file(curfile->GetName());
 257         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),"");
 258         else file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 259         TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
 260         if(!fxsec) {
 261             fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 262             if(fxsec) {
 263                 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 264                 if(key) {
 265                       TList *list = dynamic_cast<TList*>(key->ReadObj());
 266                       if(list) {
 267                           xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 268                           ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 269                       }
 270                 }
 271                 fxsec->Close();
 272             }
 273         } else {
 274             TTree *xtree = (TTree*)fxsec->Get("Xsection");
 275             if(xtree) {
 276                 UInt_t _ftrials  = 0;
 277                 Double_t _xsection  = 0;
 278                 xtree->SetBranchAddress("xsection",&_xsection);
 279                 xtree->SetBranchAddress("ntrials",&_ftrials);
 280                 xtree->GetEntry(0);
 281                 ftrials = _ftrials;
 282                 xsection = _xsection;
 283             }
 284             fxsec->Close();
 285         }
 286         fh1Xsec->Fill("<#sigma>",xsection);
 287         Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 288         if(ftrials >= nEntries && nEntries > 0.) fAvgTrials = ftrials/nEntries;
 289         fh1Trials->Fill("#sum{ntrials}",ftrials);
 290     }
 291     return kTRUE;
 292 }
 293 //_____________________________________________________________________________
 294 Bool_t AliAnalysisTaskJetMatching::Run()
 295 {
 296     // execute once for each event
 297     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 298     if(!(InputEvent() && fSourceJets && fTargetJets && IsEventSelected())) return kFALSE;
 299     if(fh1AvgTrials) fh1AvgTrials->Fill("#sum{avg ntrials}", fAvgTrials);
 300     // step one: do geometric matching
 301     switch (fMatchingScheme) {
 302         // FIXME having separate dedicated functions is not necessary and historical
 303         // cluttered code - should be cleaned up and merged into one
 304         case kGeoEtaPhi : {
 305             DoGeometricMatchingEtaPhi();
 306             // break if no jet was matched
 307             if(!fMatchedJetContainer[1][0]) return kTRUE;
 308         } break;
 309         case kGeoR : {
 310             DoGeometricMatchingR();
 311             // break if no jet was matched
 312             if(!fMatchedJetContainer[1][0]) return kTRUE;
 313         } break;
 314         case kDiJet : {
 315             DoDiJetMatching();
 316         } break;
 317         default : break;
 318     }
 319     // optional step two: get a bijection (avoid duplicate matches)
 320     if(fGetBijection)           GetBijection();
 321     // optional step three: match constituents within matched jets
 322     if(fMatchConstituents)      DoConstituentMatching();
 323     // stream data to output
 324     PostMatchedJets();
 325     FillMatchedJetHistograms();
 326     if(fDebug > 0) PrintInfo();
 327     PostData(1, fOutputList);
 328     return kTRUE;
 329 }
 330 //_____________________________________________________________________________
 331 void AliAnalysisTaskJetMatching::DoGeometricMatchingEtaPhi()
 332 {
 333     // do geometric matching based on eta phi distance between jet centers
 334     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 335     fNoMatchedJets = 0; // reset the matched jet counter
 336     Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
 337     for(Int_t i(0); i < iSource; i++) {
 338         AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
 339         if(!PassesCuts(sourceJet, 0)) continue;
 340         for(Int_t j(0); j < iTarget; j++) {
 341             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
 342             if(!PassesCuts(targetJet, 1)) continue;
 343             if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
 344             if((TMath::Abs(sourceJet->Eta() - targetJet->Eta()) < fMatchEta )) {
 345                 Double_t sourcePhi(sourceJet->Phi()), targetPhi(targetJet->Phi());
 346                 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
 347                 if(TMath::Abs(sourcePhi-targetPhi) > TMath::Abs(sourcePhi-targetPhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
 348                 if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {       // accept the jets as matching
 349                 Bool_t isBestMatch(kTRUE);
 350                     if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
 351                         if(fDebug > 0) printf(" > Entering first bijection test \n");
 352                         for(Int_t k(i); k < iSource; k++) {
 353                             AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
 354                             if(PassesCuts(candidateSourceJet, 0)) continue;
 355                             if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
 356                             if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
 357                                 isBestMatch = kFALSE;
 358                                 break;
 359                             }
 360                         }
 361                         if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
 362                     }
 363                     if(isBestMatch) {
 364                         fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
 365                         fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
 366                         fNoMatchedJets++;
 367                     }
 368                     if(fNoMatchedJets > 199) {   // maximum to keep the cache at reasonable size
 369                         AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
 370                         return;
 371                     }
 372                 }
 373             }
 374         }
 375     }
 376 }
 377 //_____________________________________________________________________________
 378 void AliAnalysisTaskJetMatching::DoGeometricMatchingR()
 379 {
 380     // do geometric matching based on shortest path between jet centers
 381     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 382     fNoMatchedJets = 0; // reset the matched jet counter
 383     Int_t iSource(fSourceJets->GetEntriesFast()), iTarget(fTargetJets->GetEntriesFast());
 384     for(Int_t i(0); i < iSource; i++) {
 385         AliEmcalJet* sourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(i)));
 386         if(!PassesCuts(sourceJet, 0)) continue;
 387         for(Int_t j(0); j < iTarget; j++) {
 388             AliEmcalJet* targetJet(static_cast<AliEmcalJet*>(fTargetJets->At(j)));
 389             if(!PassesCuts(targetJet, 1)) continue;
 390             else if (fUseEmcalBaseJetCuts && !AcceptJet(targetJet, 1)) continue;
 391             if(GetR(sourceJet, targetJet) <= fMatchR) {
 392                 Bool_t isBestMatch(kTRUE);
 393                 if(fGetBijection) { // match has been found, for bijection only store it if there's no better match
 394                     if(fDebug > 0) printf(" > Entering first bijection test \n");
 395                     for(Int_t k(i); k < iSource; k++) {
 396                         AliEmcalJet* candidateSourceJet(static_cast<AliEmcalJet*>(fSourceJets->At(k)));
 397                         if(!PassesCuts(candidateSourceJet, 0)) continue;
 398                         if(fDebug > 0) printf("source distance %.2f \t candidate distance %.2f \n", GetR(sourceJet, targetJet),GetR(candidateSourceJet, targetJet));
 399                         if(GetR(sourceJet, targetJet) > GetR(candidateSourceJet, targetJet)) {
 400                             isBestMatch = kFALSE;
 401                             break;
 402                         }
 403                     }
 404                     if(fDebug > 0) (isBestMatch) ? printf(" kept source \n ") : printf(" we can do better (rejected source) \n");
 405                 }
 406                 if(isBestMatch) {
 407                     fMatchedJetContainer[fNoMatchedJets][0] = sourceJet;
 408                     fMatchedJetContainer[fNoMatchedJets][1] = targetJet;
 409                     fNoMatchedJets++;
 410                 }
 411                 if(fNoMatchedJets > 99) {
 412                     AliError(Form("%s: Found too many matched jets (> 100). Adjust matching criteria !", GetName()));
 413                     return;
 414                 }
 415             }
 416         }
 417     }
 418 }
 419 //_____________________________________________________________________________
 420 void AliAnalysisTaskJetMatching::DoDiJetMatching()
 421 {
 422     // match dijets. this is in a sense a 'special' mode of the task as both jet supplied jet arrays
 423     // (target and source) are the same jet collection, matching will be performed on distribution in
 424     // azimuth
 425     // no ouptut array is produced, dedicated dijet histo's are filled in this function instead
 426     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 427     // get the leading jet in a given acceptance (TPC is default)
 428     Int_t leadingJetIndex(-1), subLeadingJetIndex(-1);
 429     // retrieve the leading jet, leadingJetIndex points to the leading jet
 430     AliEmcalJet* leadingJet(GetLeadingJet(fSourceJets, leadingJetIndex));
 431     if(!leadingJet) return;
 432     // fill phi and eta of leading jet (should always be in selected acceptance)
 433     fHistDiJetLeadingJet->Fill(leadingJet->Phi(), leadingJet->Eta(), leadingJet->Pt());
 434     Double_t sourcePhi(leadingJet->Phi()), targetPhi(-1);
 435     // get the sub-leading jet - faster when leading jet is also provided
 436     AliEmcalJet* subLeadingJet(GetSubLeadingJet(fSourceJets, leadingJetIndex, subLeadingJetIndex));
 437     if(!subLeadingJet) return;
 438     else { // check if the sub-leading jet is actually the away-side jet
 439         targetPhi = subLeadingJet->Phi() + TMath::Pi();
 440         // rotate jets to common phase
 441         if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi+TMath::TwoPi())) sourcePhi+=TMath::TwoPi();
 442         if(TMath::Abs(sourcePhi) > TMath::Abs(sourcePhi-TMath::TwoPi())) sourcePhi-=TMath::TwoPi();
 443         if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi+TMath::TwoPi())) targetPhi+=TMath::TwoPi();
 444         if(TMath::Abs(targetPhi) > TMath::Abs(targetPhi-TMath::TwoPi())) targetPhi-=TMath::TwoPi();
 445         if(TMath::Abs(sourcePhi-targetPhi) < fMatchPhi) {
 446             fHistDiJet->Fill(subLeadingJet->Phi(), subLeadingJet->Eta(), leadingJet->Pt());
 447             fHistDiJetDPhi->Fill(sourcePhi-targetPhi, leadingJet->Pt());
 448             fHistDiJetDPt->Fill(leadingJet->Pt() - subLeadingJet->Pt(), leadingJet->Pt());
 449         }
 450     }
 451 }
 452 //_____________________________________________________________________________
 453 void AliAnalysisTaskJetMatching::DoConstituentMatching()
 454 {
 455     // match constituents within matched jets
 456     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 457     if(!fTracks) {
 458         AliFatal(Form("%s Fatal error! To do deep matching, supply jet constituents ! \n", GetName()));
 459         return; // coverity ...
 460     }
 461     for(Int_t i(0); i < fNoMatchedJets; i++) {
 462         AliEmcalJet* sourceJet = fMatchedJetContainer[i][0];
 463         AliEmcalJet* targetJet = fMatchedJetContainer[i][1];
 464         if(sourceJet && targetJet) {    // duplicate check: slot migth be NULL
 465             Double_t targetPt(0);
 466             Int_t iSJ(sourceJet->GetNumberOfTracks());
 467             Int_t iTJ(targetJet->GetNumberOfTracks());
 468             Int_t overlap(0), alreadyFound(0);
 469             for(Int_t j(0); j < iSJ; j++) {
 470                 alreadyFound = 0;
 471                 Int_t idSource((Int_t)sourceJet->TrackAt(j));
 472                 for(Int_t k(0); k < iTJ; k++) { // compare all target tracks to the source track
 473                     if(idSource == targetJet->TrackAt(k) && alreadyFound == 0) {
 474                         overlap++;
 475                         alreadyFound++; // avoid possible duplicate matching
 476                         AliVParticle* vp(static_cast<AliVParticle*>(targetJet->TrackAt(k, fTracks)));
 477                         if(vp) targetPt += vp->Pt();
 478                         continue;
 479                     }
 480                 }
 481             }
 482             if((float)overlap/(float)iSJ < fMinFracRecoveredConstituents || targetPt/sourceJet->Pt() < fMinFracRecoveredConstituentPt) {
 483                 if(fDebug > 0) printf("  \n > Purging jet, recovered constituents ratio  %i / %i = %.2f <  or pt ratio %.2f / %.2f = %.2f < %.2f",
 484                         overlap, iSJ, (float)overlap/(float)iSJ, targetPt, sourceJet->Pt(), targetPt/sourceJet->Pt(), fMinFracRecoveredConstituentPt);
 485                 fMatchedJetContainer[i][0] = 0x0;
 486                 fMatchedJetContainer[i][1] = 0x0;
 487                 continue;
 488             }
 489             if(sourceJet->Pt() > 0) {
 490                 Double_t sourceRho(0), targetRho(0);
 491                 if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(sourceJet->Phi(), fSourceRadius)*sourceJet->Area();
 492                 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(targetJet->Phi(), fTargetRadius)*targetJet->Area();
 493                 fProfFracPtMatched->Fill(sourceJet->Pt()-sourceRho, (targetPt-targetRho) / (sourceJet->Pt()-sourceRho));
 494                 fProfFracPtJets->Fill(sourceJet->Pt()-sourceRho, (targetJet->Pt()-targetRho) / (sourceJet->Pt()-sourceRho));
 495                 fProfFracNoMatched->Fill(sourceJet->Pt()-sourceRho, (double)overlap / (double)sourceJet->GetNumberOfTracks());
 496                 fProfFracNoJets->Fill(sourceJet->Pt()-sourceRho, (double)targetJet->GetNumberOfTracks() / (double)sourceJet->GetNumberOfTracks());
 497             }
 498             if(fDebug > 0) {
 499                 printf("\n > Jet A: %i const\t", iSJ);
 500                 printf(" > Jet B %i const\t", iTJ);
 501                 printf(" -> OVERLAP: %i tracks <- \n", overlap);
 502             }
 503         }
 504     }
 505 }
 506 //_____________________________________________________________________________
 507 void AliAnalysisTaskJetMatching::GetBijection()
 508 {
 509     // bijection of source and matched jets, based on closest distance between jets
 510     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 511     for(Int_t i(0); i < fNoMatchedJets; i++) {
 512         for(Int_t j(i+1); j < fNoMatchedJets; j++) {
 513             if(fMatchedJetContainer[i][0] == fMatchedJetContainer[j][0]) {
 514                 // found source with two targets, now see which target is closer to the source
 515                 if(!(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1] && fMatchedJetContainer[j][0] && fMatchedJetContainer[j][1] )) continue;
 516                 Double_t rA(GetR(fMatchedJetContainer[i][0], fMatchedJetContainer[i][1]));      // distance between connection A = i
 517                 Double_t rB(GetR(fMatchedJetContainer[j][0], fMatchedJetContainer[j][1]));      // distance between connection B = j
 518                 if (rB > rA) {  // jet two is far away, clear it from both target and source list
 519                     fMatchedJetContainer[j][0] = 0x0;
 520                     fMatchedJetContainer[j][1] = 0x0;
 521                 } else {                // jet one is far away, clear it from both target and source list
 522                     fMatchedJetContainer[i][0] = fMatchedJetContainer[j][0];    // switch to avoid breaking loop
 523                     fMatchedJetContainer[i][1] = fMatchedJetContainer[j][1];
 524                     fMatchedJetContainer[j][0] = 0x0;                           // clear duplicate jet from cache
 525                     fMatchedJetContainer[j][1] = 0x0;
 526                 }
 527                 if(fDebug > 0) printf(" found duplicate jet, chose %.2f over %.2f \n" , (rB > rA) ? rA : rB, (rB > rA) ? rB : rA);
 528             }
 529         }
 530     }
 531 }
 532 //_____________________________________________________________________________
 533 void AliAnalysisTaskJetMatching::PostMatchedJets()
 534 {
 535     // post matched jets
 536     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 537     fMatchedJets->Delete();     // should be a NULL operation, but added just in case
 538     for(Int_t i(0), p(0); i < fNoMatchedJets; i++) {
 539         if(fMatchedJetContainer[i][1]) {        // duplicate jets will have NULL value here and are skipped
 540             new((*fMatchedJets)[p]) AliEmcalJet(*fMatchedJetContainer[i][1]);
 541             p++;
 542         }
 543     }
 544 }
 545 //_____________________________________________________________________________
 546 void AliAnalysisTaskJetMatching::FillAnalysisSummaryHistogram() const
 547 {
 548     // fill the analysis summary histrogram, saves all relevant analysis settigns
 549     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 550     fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fUseScaledRho");
 551     fHistAnalysisSummary->SetBinContent(1, (int)fUseScaledRho);
 552     fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fMatchingScheme");
 553     fHistAnalysisSummary->SetBinContent(2, (int)fMatchingScheme);
 554     fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fSourceBKG");
 555     fHistAnalysisSummary->SetBinContent(3, (int)fSourceBKG);
 556     fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fTargetBKG");
 557     fHistAnalysisSummary->SetBinContent(4, (int)fTargetBKG);
 558     fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fMatchEta");
 559     fHistAnalysisSummary->SetBinContent(5, fMatchEta);
 560     fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fMatchPhi");
 561     fHistAnalysisSummary->SetBinContent(6, fMatchPhi);
 562     fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fMatchR");
 563     fHistAnalysisSummary->SetBinContent(7, fMatchR);
 564     fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "iter");
 565     fHistAnalysisSummary->SetBinContent(8, 1.);
 566 }
 567 //_____________________________________________________________________________
 568 void AliAnalysisTaskJetMatching::FillMatchedJetHistograms()
 569 {
 570     // fill matched jet histograms
 571     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 572     for(Int_t i(0); i < fSourceJets->GetEntriesFast(); i++) {
 573         AliEmcalJet* source = static_cast<AliEmcalJet*>(fSourceJets->At(i));
 574         if(!source) continue;
 575         else if(fUseEmcalBaseJetCuts && !AcceptJet(source, 0)) continue;
 576         Double_t sourceRho(0), targetRho(0);
 577         if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(source->Phi(), fSourceRadius)*source->Area();
 578         fHistSourceJetPt->Fill(source->Pt()-sourceRho);
 579         fHistNoConstSourceJet->Fill(source->GetNumberOfConstituents());
 580         for(Int_t j(0); j < fTargetJets->GetEntriesFast(); j++) {
 581             AliEmcalJet* target = static_cast<AliEmcalJet*>(fTargetJets->At(j));
 582             if(target) {
 583                 if(fUseEmcalBaseJetCuts && !AcceptJet(target, 1)) continue;
 584                 if(fTargetRho) targetRho = fTargetRho->GetLocalVal(target->Phi(), fTargetRadius)*target->Area();
 585                 fProfQA->Fill(0.5, TMath::Abs((source->Pt()-sourceRho)-(target->Pt()-targetRho)));
 586                 fProfQA->Fill(1.5, TMath::Abs(source->Eta()-target->Eta()));
 587                 fProfQA->Fill(2.5, TMath::Abs(source->Phi()-target->Phi()));
 588                 fHistUnsortedCorrelation->Fill(PhaseShift(source->Phi()-target->Phi(), 2));
 589                 if(j==0) {
 590                     fHistTargetJetPt->Fill(target->Pt()-targetRho);
 591                     fHistNoConstTargetJet->Fill(target->GetNumberOfConstituents());
 592                 }
 593             }
 594         }
 595     }
 596     for(Int_t i(0); i < fNoMatchedJets; i++) {
 597         if(fMatchedJetContainer[i][0] && fMatchedJetContainer[i][1]) {
 598             Double_t sourceRho(0), targetRho(0);
 599             if(fSourceRho) sourceRho = fSourceRho->GetLocalVal(fMatchedJetContainer[i][0]->Phi(), fSourceRadius)*fMatchedJetContainer[i][0]->Area();
 600             if(fTargetRho) targetRho = fTargetRho->GetLocalVal(fMatchedJetContainer[i][1]->Phi(), fTargetRadius)*fMatchedJetContainer[i][1]->Area();
 601             fHistMatchedCorrelation->Fill(PhaseShift(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi(), 2));
 602             fHistMatchedSourceJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho);
 603             fHistMatchedJetPt->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
 604             fHistNoConstMatchJet->Fill(fMatchedJetContainer[i][1]->Pt()-targetRho);
 605             fProfQAMatched->Fill(0.5, TMath::Abs((fMatchedJetContainer[i][0]->Pt()-sourceRho)-(fMatchedJetContainer[i][1]->Pt()-targetRho)));
 606             fProfQAMatched->Fill(1.5, TMath::Abs(fMatchedJetContainer[i][0]->Eta()-fMatchedJetContainer[i][1]->Eta()));
 607             fProfQAMatched->Fill(2.5, TMath::Abs(fMatchedJetContainer[i][0]->Phi()-fMatchedJetContainer[i][1]->Phi()));
 608
 609             fHistSourceMatchedJetPt->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, fMatchedJetContainer[i][1]->Pt()-targetRho);
 610             if(fDoDetectorResponse) {
 611                 fProfFracPtJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (fMatchedJetContainer[i][1]->Pt()-targetRho) / (fMatchedJetContainer[i][0]->Pt()-sourceRho));
 612                 fProfFracNoJets->Fill(fMatchedJetContainer[i][0]->Pt()-sourceRho, (double)fMatchedJetContainer[i][1]->GetNumberOfTracks() / (double)fMatchedJetContainer[i][0]->GetNumberOfTracks());
 613                 fHistDetectorResponseProb->Fill((fMatchedJetContainer[i][1]->Pt()-fMatchedJetContainer[i][0]->Pt())/fMatchedJetContainer[i][0]->Pt(), fMatchedJetContainer[i][0]->Pt());
 614             }
 615         }
 616     }
 617 }
 618 //_____________________________________________________________________________
 619 AliEmcalJet* AliAnalysisTaskJetMatching::GetLeadingJet(TClonesArray* source, Int_t &leadingJetIndex, Double_t etaMin, Double_t etaMax)
 620 {
 621     // return the leading jet within giiven acceptance
 622     Int_t iJets(source->GetEntriesFast());
 623     Double_t pt(0);
 624     AliEmcalJet* leadingJet(0x0);
 625     for(Int_t i(0); i < iJets; i++) {
 626         AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
 627         if(jet->Eta() < etaMin || jet->Eta() > etaMax) continue;
 628         if(jet->Pt() > pt) {
 629            leadingJet = jet;
 630            pt = leadingJet->Pt();
 631            leadingJetIndex = i;
 632         }
 633     }
 634     return leadingJet;
 635 }
 636 //_____________________________________________________________________________
 637 AliEmcalJet* AliAnalysisTaskJetMatching::GetSubLeadingJet(TClonesArray* source, Int_t leadingJetIndex, Int_t &subLeadingJetIndex)
 638 {
 639     // return the sub-leading jet within given acceptance
 640     // same as GetLeadingJet() but skips the leading jet (so returned jet is
 641     // sub-leading by design)
 642     // there is no eta requirement on the location of the sub-leading jet
 643     Int_t iJets(source->GetEntriesFast());
 644     // if the leading jet isn't given, retrieve it
 645     if(leadingJetIndex < 0) GetLeadingJet(source, leadingJetIndex);
 646     AliEmcalJet *leadingJet(0x0), *prevLeadingJet(static_cast<AliEmcalJet*>(source->At(leadingJetIndex)));
 647     if(!prevLeadingJet) return 0x0;
 648     Double_t pt(0), leadingPt(prevLeadingJet->Pt());
 649     for(Int_t i(0); i < iJets; i++) {
 650         if(i == leadingJetIndex) continue;
 651         AliEmcalJet* jet = static_cast<AliEmcalJet*>(source->At(i));
 652         // check if jet is actually sub-leading
 653         if(jet->Pt() > pt && jet->Pt() <= leadingPt) {
 654            leadingJet = jet;
 655            pt = leadingJet->Pt();
 656            subLeadingJetIndex = i;
 657         }
 658     }
 659     return leadingJet;
 660 }
 661 //_____________________________________________________________________________
 662 void AliAnalysisTaskJetMatching::PrintInfo() const
 663 {
 664     // print some info
 665     printf("\n > No. of source jets from %s \n \t %i \n ", fSourceJetsName.Data(), fSourceJets->GetEntriesFast());
 666     printf(" > No. of target jets from %s \n \t %i \n ", fTargetJetsName.Data(), fTargetJets->GetEntriesFast());
 667     printf(" > No. of matched jets from %s \n \t %i \n ", fMatchedJetsName.Data(), fMatchedJets->GetEntriesFast());
 668     if(fDebug > 3) InputEvent()->GetList()->ls();
 669 }
 670 //_____________________________________________________________________________
 671 void AliAnalysisTaskJetMatching::Terminate(Option_t *)
 672 {
 673     // terminate
 674 }
 675 //_____________________________________________________________________________