PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx

   1 // $Id$
   2 //
   3 // Emcal jet response matrix maker task.
   4 //
   5 // Author: S. Aiola
   6
   7 #include "AliJetResponseMaker.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TH3F.h>
  13 #include <TLorentzVector.h>
  14 #include <TSystem.h>
  15 #include <TFile.h>
  16 #include <TKey.h>
  17 #include <TProfile.h>
  18
  19 #include "AliAnalysisManager.h"
  20 #include "AliVCluster.h"
  21 #include "AliVTrack.h"
  22 #include "AliEmcalJet.h"
  23 #include "AliGenPythiaEventHeader.h"
  24 #include "AliMCEvent.h"
  25 #include "AliLog.h"
  26 #include "AliRhoParameter.h"
  27 #include "AliNamedArrayI.h"
  28
  29 ClassImp(AliJetResponseMaker)
  30
  31 //________________________________________________________________________
  32 AliJetResponseMaker::AliJetResponseMaker() :
  33   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  34   fTracks2Name(""),
  35   fCalo2Name(""),
  36   fJets2Name(""),
  37   fRho2Name(""),
  38   fPtBiasJet2Track(0),
  39   fPtBiasJet2Clus(0),
  40   fAreCollections1MC(kFALSE),
  41   fAreCollections2MC(kTRUE),
  42   fMatching(kNoMatching),
  43   fMatchingPar1(0),
  44   fMatchingPar2(0),
  45   fJet2MinEta(-999),
  46   fJet2MaxEta(-999),
  47   fJet2MinPhi(-999),
  48   fJet2MaxPhi(-999),
  49   fSelectPtHardBin(-999),
  50   fPythiaHeader(0),
  51   fPtHardBin(0),
  52   fNTrials(0),
  53   fTracks2(0),
  54   fCaloClusters2(0),
  55   fJets2(0),
  56   fRho2(0),
  57   fRho2Val(0),
  58   fTracks2Map(0),
  59   fHistTrialsAfterSel(0),
  60   fHistEventsAfterSel(0),
  61   fHistTrials(0),
  62   fHistXsection(0),
  63   fHistEvents(0),
  64   fHistJets1PhiEta(0),
  65   fHistJets1PtArea(0),
  66   fHistJets1CorrPtArea(0),
  67   fHistLeadingJets1PtArea(0),
  68   fHistLeadingJets1CorrPtArea(0),
  69   fHistJets1NEFvsPt(0),
  70   fHistJets1CEFvsCEFPt(0),
  71   fHistJets1ZvsPt(0),
  72   fHistJets2PhiEta(0),
  73   fHistJets2PtArea(0),
  74   fHistJets2CorrPtArea(0),
  75   fHistLeadingJets2PtArea(0),
  76   fHistLeadingJets2CorrPtArea(0),
  77   fHistJets2PhiEtaAcceptance(0),
  78   fHistJets2PtAreaAcceptance(0),
  79   fHistJets2CorrPtAreaAcceptance(0),
  80   fHistLeadingJets2PtAreaAcceptance(0),
  81   fHistLeadingJets2CorrPtAreaAcceptance(0),
  82   fHistJets2NEFvsPt(0),
  83   fHistJets2CEFvsCEFPt(0),
  84   fHistJets2ZvsPt(0),
  85   fHistCommonEnergy1vsJet1Pt(0),
  86   fHistCommonEnergy2vsJet2Pt(0),
  87   fHistDistancevsJet1Pt(0),
  88   fHistDistancevsJet2Pt(0),
  89   fHistDistancevsCommonEnergy1(0),
  90   fHistDistancevsCommonEnergy2(0),
  91   fHistJet2PtOverJet1PtvsJet2Pt(0),
  92   fHistJet1PtOverJet2PtvsJet1Pt(0),
  93   fHistDeltaEtaPhi(0),
  94   fHistDeltaPtvsJet1Pt(0),
  95   fHistDeltaPtvsJet2Pt(0),
  96   fHistDeltaPtvsMatchingLevel(0),
  97   fHistDeltaCorrPtvsJet1Pt(0),
  98   fHistDeltaCorrPtvsJet2Pt(0),
  99   fHistDeltaCorrPtvsMatchingLevel(0),
 100   fHistNonMatchedJets1PtArea(0),
 101   fHistNonMatchedJets2PtArea(0),
 102   fHistNonMatchedJets1CorrPtArea(0),
 103   fHistNonMatchedJets2CorrPtArea(0),
 104   fHistJet1PtvsJet2Pt(0),
 105   fHistJet1CorrPtvsJet2CorrPt(0),
 106   fHistMissedJets2PtArea(0)
 107 {
 108   // Default constructor.
 109
 110   SetMakeGeneralHistograms(kTRUE);
 111 }
 112
 113 //________________________________________________________________________
 114 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
 115   AliAnalysisTaskEmcalJet(name, kTRUE),
 116   fTracks2Name("MCParticles"),
 117   fCalo2Name(""),
 118   fJets2Name("MCJets"),
 119   fRho2Name(""),
 120   fPtBiasJet2Track(0),
 121   fPtBiasJet2Clus(0),
 122   fAreCollections1MC(kFALSE),
 123   fAreCollections2MC(kTRUE),
 124   fMatching(kNoMatching),
 125   fMatchingPar1(0),
 126   fMatchingPar2(0),
 127   fJet2MinEta(-999),
 128   fJet2MaxEta(-999),
 129   fJet2MinPhi(-999),
 130   fJet2MaxPhi(-999),
 131   fSelectPtHardBin(-999),
 132   fPythiaHeader(0),
 133   fPtHardBin(0),
 134   fNTrials(0),
 135   fTracks2(0),
 136   fCaloClusters2(0),
 137   fJets2(0),
 138   fRho2(0),
 139   fRho2Val(0),
 140   fTracks2Map(0),
 141   fHistTrialsAfterSel(0),
 142   fHistEventsAfterSel(0),
 143   fHistTrials(0),
 144   fHistXsection(0),
 145   fHistEvents(0),
 146   fHistJets1PhiEta(0),
 147   fHistJets1PtArea(0),
 148   fHistJets1CorrPtArea(0),
 149   fHistLeadingJets1PtArea(0),
 150   fHistLeadingJets1CorrPtArea(0),
 151   fHistJets1NEFvsPt(0),
 152   fHistJets1CEFvsCEFPt(0),
 153   fHistJets1ZvsPt(0),
 154   fHistJets2PhiEta(0),
 155   fHistJets2PtArea(0),
 156   fHistJets2CorrPtArea(0),
 157   fHistLeadingJets2PtArea(0),
 158   fHistLeadingJets2CorrPtArea(0),
 159   fHistJets2PhiEtaAcceptance(0),
 160   fHistJets2PtAreaAcceptance(0),
 161   fHistJets2CorrPtAreaAcceptance(0),
 162   fHistLeadingJets2PtAreaAcceptance(0),
 163   fHistLeadingJets2CorrPtAreaAcceptance(0),
 164   fHistJets2NEFvsPt(0),
 165   fHistJets2CEFvsCEFPt(0),
 166   fHistJets2ZvsPt(0),
 167   fHistCommonEnergy1vsJet1Pt(0),
 168   fHistCommonEnergy2vsJet2Pt(0),
 169   fHistDistancevsJet1Pt(0),
 170   fHistDistancevsJet2Pt(0),
 171   fHistDistancevsCommonEnergy1(0),
 172   fHistDistancevsCommonEnergy2(0),
 173   fHistJet2PtOverJet1PtvsJet2Pt(0),
 174   fHistJet1PtOverJet2PtvsJet1Pt(0),
 175   fHistDeltaEtaPhi(0),
 176   fHistDeltaPtvsJet1Pt(0),
 177   fHistDeltaPtvsJet2Pt(0),
 178   fHistDeltaPtvsMatchingLevel(0),
 179   fHistDeltaCorrPtvsJet1Pt(0),
 180   fHistDeltaCorrPtvsJet2Pt(0),
 181   fHistDeltaCorrPtvsMatchingLevel(0),
 182   fHistNonMatchedJets1PtArea(0),
 183   fHistNonMatchedJets2PtArea(0),
 184   fHistNonMatchedJets1CorrPtArea(0),
 185   fHistNonMatchedJets2CorrPtArea(0),
 186   fHistJet1PtvsJet2Pt(0),
 187   fHistJet1CorrPtvsJet2CorrPt(0),
 188   fHistMissedJets2PtArea(0)
 189 {
 190   // Standard constructor.
 191
 192   SetMakeGeneralHistograms(kTRUE);
 193 }
 194
 195 //________________________________________________________________________
 196 AliJetResponseMaker::~AliJetResponseMaker()
 197 {
 198   // Destructor
 199 }
 200
 201 //________________________________________________________________________
 202 Bool_t AliJetResponseMaker::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
 203 {
 204   //
 205   // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 206   // Get the pt hard bin from the file path
 207   // This is to called in Notify and should provide the path to the AOD/ESD file
 208   // (Partially copied from AliAnalysisHelperJetTasks)
 209
 210   TString file(currFile);
 211   fXsec = 0;
 212   fTrials = 1;
 213
 214   if(file.Contains(".zip#")){
 215     Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
 216     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 217     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
 218     file.Replace(pos+1,pos2-pos1,"");
 219   }
 220   else {
 221     // not an archive take the basename....
 222     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 223   }
 224   Printf("%s",file.Data());
 225
 226   // Get the pt hard bin
 227   TString strPthard(file);
 228   strPthard.Remove(strPthard.Last('/'));
 229   strPthard.Remove(strPthard.Last('/'));
 230   strPthard.Remove(0,strPthard.Last('/')+1);
 231   if (strPthard.IsDec())
 232     pthard = strPthard.Atoi();
 233   else
 234     AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
 235
 236   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 237   if(!fxsec){
 238     // next trial fetch the histgram file
 239     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 240     if(!fxsec){
 241         // not a severe condition but inciate that we have no information
 242       return kFALSE;
 243     }
 244     else{
 245       // find the tlist we want to be independtent of the name so use the Tkey
 246       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 247       if(!key){
 248         fxsec->Close();
 249         return kFALSE;
 250       }
 251       TList *list = dynamic_cast<TList*>(key->ReadObj());
 252       if(!list){
 253         fxsec->Close();
 254         return kFALSE;
 255       }
 256       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 257       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 258       fxsec->Close();
 259     }
 260   } // no tree pyxsec.root
 261   else {
 262     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 263     if(!xtree){
 264       fxsec->Close();
 265       return kFALSE;
 266     }
 267     UInt_t   ntrials  = 0;
 268     Double_t  xsection  = 0;
 269     xtree->SetBranchAddress("xsection",&xsection);
 270     xtree->SetBranchAddress("ntrials",&ntrials);
 271     xtree->GetEntry(0);
 272     fTrials = ntrials;
 273     fXsec = xsection;
 274     fxsec->Close();
 275   }
 276   return kTRUE;
 277 }
 278
 279 //________________________________________________________________________
 280 Bool_t AliJetResponseMaker::UserNotify()
 281 {
 282   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 283   if (!tree) {
 284     AliError(Form("%s - UserNotify: No current tree!",GetName()));
 285     return kFALSE;
 286   }
 287
 288   Float_t xsection = 0;
 289   Float_t trials   = 0;
 290   Int_t   pthard   = 0;
 291
 292   TFile *curfile = tree->GetCurrentFile();
 293   if (!curfile) {
 294     AliError(Form("%s - UserNotify: No current file!",GetName()));
 295     return kFALSE;
 296   }
 297
 298   PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
 299
 300   fHistTrials->SetBinContent(pthard + 1, fHistTrials->GetBinContent(pthard + 1) + trials);
 301   fHistXsection->SetBinContent(pthard + 1, fHistXsection->GetBinContent(pthard + 1) + xsection);
 302   fHistEvents->SetBinContent(pthard + 1, fHistEvents->GetBinContent(pthard + 1) + 1);
 303
 304   return kTRUE;
 305 }
 306
 307 //________________________________________________________________________
 308 void AliJetResponseMaker::UserCreateOutputObjects()
 309 {
 310   // Create user objects.
 311
 312   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 313
 314   fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
 315   fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
 316   fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
 317   fOutput->Add(fHistTrialsAfterSel);
 318
 319   fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
 320   fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
 321   fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
 322   fOutput->Add(fHistEventsAfterSel);
 323
 324   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
 325   fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 326   fHistTrials->GetYaxis()->SetTitle("trials");
 327   fOutput->Add(fHistTrials);
 328
 329   fHistXsection = new TH1F("fHistXsection", "fHistXsection", 11, 0, 11);
 330   fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
 331   fHistXsection->GetYaxis()->SetTitle("xsection");
 332   fOutput->Add(fHistXsection);
 333
 334   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 335   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 336   fHistEvents->GetYaxis()->SetTitle("total events");
 337   fOutput->Add(fHistEvents);
 338
 339   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 340   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 341
 342   for (Int_t i = 1; i < 12; i++) {
 343     fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 344     fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 345
 346     fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 347     fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 348     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 349   }
 350
 351   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 352   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
 353   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
 354   fOutput->Add(fHistJets1PhiEta);
 355
 356   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 357   fHistJets1PtArea->GetXaxis()->SetTitle("area");
 358   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 359   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
 360   fOutput->Add(fHistJets1PtArea);
 361
 362   fHistLeadingJets1PtArea = new TH2F("fHistLeadingJets1PtArea", "fHistLeadingJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 363   fHistLeadingJets1PtArea->GetXaxis()->SetTitle("area");
 364   fHistLeadingJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 365   fHistLeadingJets1PtArea->GetZaxis()->SetTitle("counts");
 366   fOutput->Add(fHistLeadingJets1PtArea);
 367
 368   if (!fRhoName.IsNull()) {
 369     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 370     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
 371     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
 372     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 373     fOutput->Add(fHistJets1CorrPtArea);
 374
 375     fHistLeadingJets1CorrPtArea = new TH2F("fHistLeadingJets1CorrPtArea", "fHistLeadingJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 376     fHistLeadingJets1CorrPtArea->GetXaxis()->SetTitle("area");
 377     fHistLeadingJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
 378     fHistLeadingJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 379     fOutput->Add(fHistLeadingJets1CorrPtArea);
 380   }
 381
 382   fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 383   fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
 384   fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 385   fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
 386   fOutput->Add(fHistJets1ZvsPt);
 387
 388   fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 389   fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
 390   fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 391   fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
 392   fOutput->Add(fHistJets1NEFvsPt);
 393
 394   fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 395   fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
 396   fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
 397   fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
 398   fOutput->Add(fHistJets1CEFvsCEFPt);
 399
 400   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 401   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
 402   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
 403   fOutput->Add(fHistJets2PhiEta);
 404
 405   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
 406   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
 407   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
 408   fOutput->Add(fHistJets2PhiEtaAcceptance);
 409
 410   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 411   fHistJets2PtArea->GetXaxis()->SetTitle("area");
 412   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 413   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
 414   fOutput->Add(fHistJets2PtArea);
 415
 416   fHistLeadingJets2PtArea = new TH2F("fHistLeadingJets2PtArea", "fHistLeadingJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 417   fHistLeadingJets2PtArea->GetXaxis()->SetTitle("area");
 418   fHistLeadingJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 419   fHistLeadingJets2PtArea->GetZaxis()->SetTitle("counts");
 420   fOutput->Add(fHistLeadingJets2PtArea);
 421
 422   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 423   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
 424   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 425   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
 426   fOutput->Add(fHistJets2PtAreaAcceptance);
 427
 428   fHistLeadingJets2PtAreaAcceptance = new TH2F("fHistLeadingJets2PtAreaAcceptance", "fHistLeadingJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 429   fHistLeadingJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
 430   fHistLeadingJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 431   fHistLeadingJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
 432   fOutput->Add(fHistLeadingJets2PtAreaAcceptance);
 433
 434   if (!fRho2Name.IsNull()) {
 435     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 436     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
 437     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 438     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 439     fOutput->Add(fHistJets2CorrPtArea);
 440
 441     fHistLeadingJets2CorrPtArea = new TH2F("fHistLeadingJets2CorrPtArea", "fHistLeadingJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 442     fHistLeadingJets2CorrPtArea->GetXaxis()->SetTitle("area");
 443     fHistLeadingJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 444     fHistLeadingJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 445     fOutput->Add(fHistLeadingJets2CorrPtArea);
 446
 447     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 448     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
 449     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 450     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
 451     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
 452
 453     fHistLeadingJets2CorrPtAreaAcceptance = new TH2F("fHistLeadingJets2CorrPtAreaAcceptance", "fHistLeadingJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 454     fHistLeadingJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
 455     fHistLeadingJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 456     fHistLeadingJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
 457     fOutput->Add(fHistLeadingJets2CorrPtAreaAcceptance);
 458   }
 459
 460   fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 461   fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
 462   fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 463   fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
 464   fOutput->Add(fHistJets2ZvsPt);
 465
 466   fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 467   fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
 468   fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 469   fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
 470   fOutput->Add(fHistJets2NEFvsPt);
 471
 472   fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 473   fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
 474   fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
 475   fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
 476   fOutput->Add(fHistJets2CEFvsCEFPt);
 477
 478   fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 479   fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
 480   fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
 481   fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
 482   fOutput->Add(fHistCommonEnergy1vsJet1Pt);
 483
 484   fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 485   fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
 486   fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 487   fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
 488   fOutput->Add(fHistCommonEnergy2vsJet2Pt);
 489
 490   fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 491   fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
 492   fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
 493   fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
 494   fOutput->Add(fHistDistancevsJet1Pt);
 495
 496   fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 497   fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
 498   fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 499   fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
 500   fOutput->Add(fHistDistancevsJet2Pt);
 501
 502   fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 503   fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
 504   fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
 505   fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
 506   fOutput->Add(fHistDistancevsCommonEnergy1);
 507
 508   fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 509   fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
 510   fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
 511   fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
 512   fOutput->Add(fHistDistancevsCommonEnergy2);
 513
 514   fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
 515   fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 516   fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
 517   fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 518   fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
 519
 520   fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
 521   fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 522   fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
 523   fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
 524   fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
 525
 526   fHistDeltaEtaPhi = new TH2F("fHistDeltaEtaPhi", "fHistDeltaEtaPhi", 200, -1, 1, 250, -1.6, 4.8);
 527   fHistDeltaEtaPhi->GetXaxis()->SetTitle("#Delta#eta");
 528   fHistDeltaEtaPhi->GetYaxis()->SetTitle("#Delta#phi");
 529   fHistDeltaEtaPhi->GetZaxis()->SetTitle("counts");
 530   fOutput->Add(fHistDeltaEtaPhi);
 531
 532   fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 533   fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 534   fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 535   fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
 536   fOutput->Add(fHistDeltaPtvsJet1Pt);
 537
 538   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 539   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 540   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 541   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 542   fOutput->Add(fHistDeltaPtvsJet2Pt);
 543
 544   fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 545   fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
 546   fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
 547   fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
 548   fOutput->Add(fHistDeltaPtvsMatchingLevel);
 549
 550   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 551     fHistDeltaCorrPtvsJet1Pt = new TH2F("fHistDeltaCorrPtvsJet1Pt", "fHistDeltaCorrPtvsJet1Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 552     fHistDeltaCorrPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 553     fHistDeltaCorrPtvsJet1Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 554     fHistDeltaCorrPtvsJet1Pt->GetZaxis()->SetTitle("counts");
 555     fOutput->Add(fHistDeltaCorrPtvsJet1Pt);
 556
 557     fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 558     fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 559     fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 560     fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 561     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
 562
 563     fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 564     fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
 565     fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 566     fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
 567     fOutput->Add(fHistDeltaCorrPtvsMatchingLevel);
 568   }
 569
 570   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 571   fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
 572   fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 573   fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
 574   fOutput->Add(fHistNonMatchedJets1PtArea);
 575
 576   fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 577   fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
 578   fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 579   fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
 580   fOutput->Add(fHistNonMatchedJets2PtArea);
 581
 582   if (!fRhoName.IsNull()) {
 583     fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 584     fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
 585     fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 586     fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 587     fOutput->Add(fHistNonMatchedJets1CorrPtArea);
 588   }
 589
 590   if (!fRho2Name.IsNull()) {
 591     fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 592     fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
 593     fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 594     fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 595     fOutput->Add(fHistNonMatchedJets2CorrPtArea);
 596   }
 597
 598   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 599   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
 600   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 601   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 602   fOutput->Add(fHistJet1PtvsJet2Pt);
 603
 604   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 605     if (fRhoName.IsNull())
 606       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 607     else if (fRho2Name.IsNull())
 608       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 609     else
 610       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 611     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 612     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
 613     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 614     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
 615   }
 616
 617   fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 618   fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
 619   fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 620   fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
 621   fOutput->Add(fHistMissedJets2PtArea);
 622
 623   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 624 }
 625
 626 //________________________________________________________________________
 627 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 628 {
 629   // Return true if jet is accepted.
 630
 631   if (jet->Pt() <= fJetPtCut)
 632     return kFALSE;
 633   if (jet->Area() <= fJetAreaCut)
 634     return kFALSE;
 635
 636   return kTRUE;
 637 }
 638
 639 //________________________________________________________________________
 640 Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
 641 {
 642   // Accept jet with a bias.
 643
 644   if (fLeadingHadronType == 0) {
 645     if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
 646   }
 647   else if (fLeadingHadronType == 1) {
 648     if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 649   }
 650   else {
 651     if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 652   }
 653
 654   return kTRUE;
 655 }
 656
 657 //________________________________________________________________________
 658 void AliJetResponseMaker::ExecOnce()
 659 {
 660   // Execute once.
 661
 662   if (!fJets2Name.IsNull() && !fJets2) {
 663     fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
 664     if (!fJets2) {
 665       AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
 666       return;
 667     }
 668     else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
 669       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
 670       fJets2 = 0;
 671       return;
 672     }
 673   }
 674
 675   if (!fTracks2Name.IsNull() && !fTracks2) {
 676     fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
 677     if (!fTracks2) {
 678       AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
 679       return;
 680     }
 681     else {
 682       TClass *cl = fTracks2->GetClass();
 683       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 684         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
 685         fTracks2 = 0;
 686         return;
 687       }
 688     }
 689
 690     if (fAreCollections2MC) {
 691       fTracks2Map = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
 692       // this is needed to map the MC labels with the indexes of the MC particle collection
 693       // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
 694       if (!fTracks2Map) {
 695         AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data()));
 696         fTracks2Map = new AliNamedArrayI("tracksMap",9999);
 697         for (Int_t i = 0; i < 9999; i++) {
 698           fTracks2Map->AddAt(i,i);
 699         }
 700       }
 701     }
 702   }
 703
 704   if (!fCalo2Name.IsNull() && !fCaloClusters2) {
 705     fCaloClusters2 =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 706     if (!fCaloClusters2) {
 707       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
 708       return;
 709     } else {
 710       TClass *cl = fCaloClusters2->GetClass();
 711       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
 712         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
 713         fCaloClusters2 = 0;
 714         return;
 715       }
 716     }
 717   }
 718
 719   if (!fRho2Name.IsNull() && !fRho2) {
 720     fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
 721     if (!fRho2) {
 722       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
 723       fInitialized = kFALSE;
 724       return;
 725     }
 726   }
 727
 728   if (fJet2MinEta == -999)
 729     fJet2MinEta = fJetMinEta - fJetRadius;
 730   if (fJet2MaxEta == -999)
 731     fJet2MaxEta = fJetMaxEta + fJetRadius;
 732   if (fJet2MinPhi == -999)
 733     fJet2MinPhi = fJetMinPhi - fJetRadius;
 734   if (fJet2MaxPhi == -999)
 735     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
 736
 737   AliAnalysisTaskEmcalJet::ExecOnce();
 738 }
 739
 740 //________________________________________________________________________
 741 Bool_t AliJetResponseMaker::IsEventSelected()
 742 {
 743   // Check if event is selected
 744
 745   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
 746     return kFALSE;
 747
 748   return AliAnalysisTaskEmcalJet::IsEventSelected();
 749 }
 750
 751 //________________________________________________________________________
 752 Bool_t AliJetResponseMaker::RetrieveEventObjects()
 753 {
 754   // Retrieve event objects.
 755
 756   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 757     return kFALSE;
 758
 759   if (fRho2)
 760     fRho2Val = fRho2->GetVal();
 761
 762   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 763   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 764
 765   if (MCEvent())
 766     fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 767
 768   if (fPythiaHeader) {
 769     Double_t pthard = fPythiaHeader->GetPtHard();
 770
 771     for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 772       if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
 773         break;
 774     }
 775
 776     fNTrials = fPythiaHeader->Trials();
 777   }
 778
 779   return kTRUE;
 780 }
 781
 782 //________________________________________________________________________
 783 Bool_t AliJetResponseMaker::Run()
 784 {
 785   // Find the closest jets
 786
 787   if (fMatching == kNoMatching)
 788     return kTRUE;
 789   else
 790     return DoJetMatching();
 791 }
 792
 793 //________________________________________________________________________
 794 Bool_t AliJetResponseMaker::DoJetMatching()
 795 {
 796   DoJetLoop(kFALSE);
 797
 798   const Int_t nJets = fJets->GetEntriesFast();
 799
 800   for (Int_t i = 0; i < nJets; i++) {
 801
 802     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
 803
 804     if (!jet1) {
 805       AliError(Form("Could not receive jet %d", i));
 806       continue;
 807     }
 808
 809     if (!AcceptJet(jet1))
 810       continue;
 811
 812     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 813       continue;
 814
 815     if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 &&
 816         jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) {    // Matched jet found
 817       jet1->SetMatchedToClosest(fMatching);
 818       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
 819     }
 820   }
 821
 822   return kTRUE;
 823 }
 824
 825 //________________________________________________________________________
 826 void AliJetResponseMaker::DoJetLoop(Bool_t order)
 827 {
 828   // Do the jet loop.
 829
 830   TClonesArray *jets1 = 0;
 831   TClonesArray *jets2 = 0;
 832
 833   if (order) {
 834     jets1 = fJets2;
 835     jets2 = fJets;
 836   }
 837   else {
 838     jets1 = fJets;
 839     jets2 = fJets2;
 840   }
 841
 842   Int_t nJets1 = jets1->GetEntriesFast();
 843   Int_t nJets2 = jets2->GetEntriesFast();
 844
 845   Bool_t jet2Reset = kFALSE;
 846
 847   for (Int_t i = 0; i < nJets1; i++) {
 848
 849     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 850
 851     if (!jet1) {
 852       AliError(Form("Could not receive jet %d", i));
 853       continue;
 854     }
 855
 856     jet1->ResetMatching();
 857
 858     if (!AcceptJet(jet1))
 859       continue;
 860
 861     if (order) {
 862      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 863         continue;
 864     }
 865     else {
 866       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 867         continue;
 868     }
 869
 870     for (Int_t j = 0; j < nJets2; j++) {
 871
 872       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
 873
 874       if (!jet2) {
 875         AliError(Form("Could not receive jet %d", j));
 876         continue;
 877       }
 878
 879       if (!jet2Reset)
 880         jet2->ResetMatching();
 881
 882       if (!AcceptJet(jet2))
 883         continue;
 884
 885       if (order) {
 886         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
 887           continue;
 888       }
 889       else {
 890         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 891           continue;
 892       }
 893
 894       SetMatchingLevel(jet1, jet2, fMatching);
 895     } // jet2 loop
 896
 897     jet2Reset = kTRUE;
 898   } // jet1 loop
 899 }
 900
 901 //________________________________________________________________________
 902 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
 903 {
 904   Double_t deta = jet2->Eta() - jet1->Eta();
 905   Double_t dphi = jet2->Phi() - jet1->Phi();
 906   d = TMath::Sqrt(deta * deta + dphi * dphi);
 907 }
 908
 909 //________________________________________________________________________
 910 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
 911 {
 912   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
 913   d1 = jet1->Pt();
 914   d2 = jet2->Pt();
 915   Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
 916
 917   for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
 918     Bool_t track2Found = kFALSE;
 919     Int_t index2 = jet2->TrackAt(iTrack2);
 920     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
 921       AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
 922       if (!track) {
 923         AliWarning(Form("Could not find track %d!", iTrack));
 924         continue;
 925       }
 926       Int_t MClabel = TMath::Abs(track->GetLabel());
 927       Int_t index = -1;
 928
 929       if (MClabel == 0) {// this is not a MC particle; remove it completely
 930         AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
 931         totalPt1 -= track->Pt();
 932         d1 -= track->Pt();
 933         continue;
 934       }
 935       else if (MClabel < fTracks2Map->GetSize()) {
 936         index = fTracks2Map->At(MClabel);
 937       }
 938
 939       if (index < 0) {
 940         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
 941         continue;
 942       }
 943
 944       if (index2 == index) { // found common particle
 945         track2Found = kTRUE;
 946         d1 -= track->Pt();
 947         AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
 948         AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
 949                         iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
 950         d2 -= MCpart->Pt();
 951         break;
 952       }
 953     }
 954     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
 955       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
 956       if (!clus) {
 957         AliWarning(Form("Could not find cluster %d!", iClus));
 958         continue;
 959       }
 960       TLorentzVector part;
 961       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
 962
 963       if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
 964         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
 965           Int_t cellId = clus->GetCellAbsId(iCell);
 966           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
 967
 968           Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
 969           Int_t index = -1;
 970
 971           if (MClabel == 0) {// this is not a MC particle; remove it completely
 972             AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
 973             totalPt1 -= part.Pt() * cellFrac;
 974             d1 -= part.Pt() * cellFrac;
 975             continue;
 976           }
 977           else if (MClabel < fTracks2Map->GetSize()) {
 978             index = fTracks2Map->At(MClabel);
 979           }
 980
 981           if (index < 0) {
 982             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
 983             continue;
 984           }
 985           if (index2 == index) { // found common particle
 986             d1 -= part.Pt() * cellFrac;
 987
 988             if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
 989               AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
 990               AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
 991                               iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
 992               d2 -= MCpart->Pt() * cellFrac;
 993             }
 994             break;
 995           }
 996         }
 997       }
 998       else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
 999         Int_t MClabel = TMath::Abs(clus->GetLabel());
1000         Int_t index = -1;
1001
1002         if (MClabel == 0) {// this is not a MC particle; remove it completely
1003           AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1004           totalPt1 -= part.Pt();
1005           d1 -= part.Pt();
1006           continue;
1007         }
1008         else if (MClabel < fTracks2Map->GetSize()) {
1009           index = fTracks2Map->At(MClabel);
1010         }
1011
1012         if (index < 0) {
1013           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1014           continue;
1015         }
1016         if (index2 == index) { // found common particle
1017           d1 -= part.Pt();
1018
1019           if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
1020             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
1021             AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1022                             iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1023
1024             d2 -= MCpart->Pt();
1025           }
1026           break;
1027         }
1028       }
1029     }
1030   }
1031   if (d1 <= 0 || totalPt1 < 1)
1032     d1 = 0;
1033   else
1034     d1 /= totalPt1;
1035
1036   if (jet2->Pt() > 0 && d2 > 0)
1037     d2 /= jet2->Pt();
1038   else
1039     d2 = 0;
1040 }
1041
1042 //________________________________________________________________________
1043 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1044 {
1045   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1046   d1 = jet1->Pt();
1047   d2 = jet2->Pt();
1048
1049   if (fTracks && fTracks2) {
1050
1051     for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1052       Int_t index2 = jet2->TrackAt(iTrack2);
1053       for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1054         Int_t index = jet1->TrackAt(iTrack);
1055         if (index2 == index) { // found common particle
1056           AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
1057           if (!part) {
1058             AliWarning(Form("Could not find track %d!", index));
1059             continue;
1060           }
1061           AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
1062           if (!part2) {
1063             AliWarning(Form("Could not find track %d!", index2));
1064             continue;
1065           }
1066
1067           d1 -= part->Pt();
1068           d2 -= part2->Pt();
1069           break;
1070         }
1071       }
1072     }
1073
1074   }
1075
1076   if (fCaloClusters && fCaloClusters2) {
1077
1078     for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1079       Int_t index2 = jet2->ClusterAt(iClus2);
1080       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1081         Int_t index = jet1->ClusterAt(iClus);
1082         AliVCluster *clus =  static_cast<AliVCluster*>(fCaloClusters->At(index));
1083         if (!clus) {
1084           AliWarning(Form("Could not find cluster %d!", index));
1085           continue;
1086         }
1087         AliVCluster *clus2 =  static_cast<AliVCluster*>(fCaloClusters2->At(index2));
1088         if (!clus2) {
1089           AliWarning(Form("Could not find cluster %d!", index2));
1090           continue;
1091         }
1092         TLorentzVector part, part2;
1093         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1094         clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1095
1096         d1 -= part.Pt();
1097         d2 -= part2.Pt();
1098         break;
1099       }
1100     }
1101
1102   }
1103
1104   if (jet1->Pt() > 0 && d1 > 0)
1105     d1 /= jet1->Pt();
1106   else
1107     d1 = 0;
1108
1109   if (jet2->Pt() > 0 && d2 > 0)
1110     d2 /= jet2->Pt();
1111   else
1112     d2 = 0;
1113 }
1114
1115 //________________________________________________________________________
1116 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
1117 {
1118   Double_t d1 = -1;
1119   Double_t d2 = -1;
1120
1121   switch (matching) {
1122   case kGeometrical:
1123     GetGeometricalMatchingLevel(jet1,jet2,d1);
1124     d2 = d1;
1125     break;
1126   case kMCLabel: // jet1 = detector level and jet2 = particle level!
1127     GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1128     break;
1129   case kSameCollections:
1130     GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1131     break;
1132   default:
1133     ;
1134   }
1135
1136   if (d1 > 0) {
1137
1138     if (d1 < jet1->ClosestJetDistance()) {
1139       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1140       jet1->SetClosestJet(jet2, d1);
1141     }
1142     else if (d1 < jet1->SecondClosestJetDistance()) {
1143       jet1->SetSecondClosestJet(jet2, d1);
1144     }
1145   }
1146
1147   if (d2 > 0) {
1148
1149     if (d2 < jet2->ClosestJetDistance()) {
1150       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1151       jet2->SetClosestJet(jet1, d2);
1152     }
1153     else if (d2 < jet2->SecondClosestJetDistance()) {
1154       jet2->SetSecondClosestJet(jet1, d2);
1155     }
1156   }
1157 }
1158
1159 //________________________________________________________________________
1160 Bool_t AliJetResponseMaker::FillHistograms()
1161 {
1162   // Fill histograms.
1163
1164   static Int_t indexes[9999] = {-1};
1165
1166   fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
1167   fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
1168
1169   GetSortedArray(indexes, fJets2, fRho2Val);
1170
1171   const Int_t nJets2 = fJets2->GetEntriesFast();
1172
1173   Int_t naccJets2 = 0;
1174   Int_t naccJets2Acceptance = 0;
1175
1176   for (Int_t i = 0; i < nJets2; i++) {
1177
1178     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(indexes[i]));
1179
1180     if (!jet2) {
1181       AliError(Form("Could not receive jet2 %d", i));
1182       continue;
1183     }
1184
1185     if (!AcceptJet(jet2))
1186       continue;
1187
1188     if (AcceptBiasJet(jet2) &&
1189         (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
1190
1191       fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
1192       fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
1193
1194       if (naccJets2Acceptance < fNLeadingJets)
1195         fHistLeadingJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
1196
1197       if (!fRho2Name.IsNull()) {
1198         fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1199         if (naccJets2Acceptance < fNLeadingJets)
1200           fHistLeadingJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1201       }
1202
1203       if (fTracks2) {
1204         for (Int_t it = 0; it < jet2->GetNumberOfTracks(); it++) {
1205           AliVParticle *track2 = jet2->TrackAt(it, fTracks2);
1206           if (track2)
1207             fHistJets2ZvsPt->Fill(track2->Pt() / jet2->Pt(), jet2->Pt());
1208         }
1209       }
1210
1211       if (fCaloClusters2) {
1212         for (Int_t ic = 0; ic < jet2->GetNumberOfClusters(); ic++) {
1213           AliVCluster *cluster2 = jet2->ClusterAt(ic, fCaloClusters2);
1214
1215           if (cluster2) {
1216             TLorentzVector nPart2;
1217             cluster2->GetMomentum(nPart2, fVertex);
1218             fHistJets2ZvsPt->Fill(nPart2.Et() / jet2->Pt(), jet2->Pt());
1219           }
1220         }
1221       }
1222
1223       fHistJets2NEFvsPt->Fill(jet2->NEF(), jet2->Pt());
1224       fHistJets2CEFvsCEFPt->Fill(1-jet2->NEF(), (1-jet2->NEF())*jet2->Pt());
1225
1226       naccJets2Acceptance++;
1227     }
1228
1229     if (!AcceptBiasJet2(jet2))
1230       continue;
1231
1232     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
1233       continue;
1234
1235     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1236     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
1237
1238     if (naccJets2 < fNLeadingJets)
1239       fHistLeadingJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1240
1241     if (!fRho2Name.IsNull()) {
1242       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1243       if (naccJets2 < fNLeadingJets)
1244         fHistLeadingJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1245     }
1246
1247     naccJets2++;
1248
1249     if (jet2->MatchedJet()) {
1250
1251       if (!AcceptBiasJet(jet2->MatchedJet()) ||
1252           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt) {
1253         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1254       }
1255       else {
1256         if (jet2->MatchedJet()->Pt() > fMaxBinPt)
1257           fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1258
1259         Double_t d1=-1, d2=-1;
1260         if (jet2->GetMatchingType() == kGeometrical) {
1261
1262           if (fAreCollections2MC && !fAreCollections1MC)
1263             GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
1264           else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
1265             GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
1266
1267           fHistDistancevsCommonEnergy1->Fill(jet2->ClosestJetDistance(), d1);
1268           fHistDistancevsCommonEnergy2->Fill(jet2->ClosestJetDistance(), d2);
1269
1270           fHistDistancevsJet1Pt->Fill(jet2->ClosestJetDistance(), jet2->MatchedJet()->Pt());
1271           fHistDistancevsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
1272
1273           fHistCommonEnergy1vsJet1Pt->Fill(d1, jet2->MatchedJet()->Pt());
1274           fHistCommonEnergy2vsJet2Pt->Fill(d2, jet2->Pt());
1275         }
1276         else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1277           GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
1278
1279           fHistDistancevsCommonEnergy1->Fill(d1, jet2->MatchedJet()->ClosestJetDistance());
1280           fHistDistancevsCommonEnergy2->Fill(d1, jet2->ClosestJetDistance());
1281
1282           fHistDistancevsJet1Pt->Fill(d1, jet2->MatchedJet()->Pt());
1283           fHistDistancevsJet2Pt->Fill(d1, jet2->Pt());
1284
1285           fHistCommonEnergy1vsJet1Pt->Fill(jet2->MatchedJet()->ClosestJetDistance(), jet2->MatchedJet()->Pt());
1286           fHistCommonEnergy2vsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
1287         }
1288
1289         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
1290         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
1291         fHistDeltaEtaPhi->Fill(deta, dphi, jet2->Pt());
1292
1293         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
1294         fHistDeltaPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
1295         fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
1296         fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
1297
1298         fHistJet2PtOverJet1PtvsJet2Pt->Fill(jet2->Pt(), jet2->Pt() / jet2->MatchedJet()->Pt());
1299         fHistJet1PtOverJet2PtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), jet2->MatchedJet()->Pt() / jet2->Pt());
1300
1301         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
1302
1303         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
1304           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
1305           fHistDeltaCorrPtvsJet1Pt->Fill(jet2->MatchedJet()->Pt(), dpt);
1306           fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
1307           fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
1308           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
1309         }
1310       }
1311     }
1312     else {
1313       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1314       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
1315
1316       if (!fRho2Name.IsNull())
1317         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
1318     }
1319   }
1320
1321   GetSortedArray(indexes, fJets, fRhoVal);
1322
1323   const Int_t nJets1 = fJets->GetEntriesFast();
1324   Int_t naccJets1 = 0;
1325
1326   for (Int_t i = 0; i < nJets1; i++) {
1327
1328     AliDebug(2,Form("Processing jet %d", i));
1329     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(indexes[i]));
1330
1331     if (!jet1) {
1332       AliError(Form("Could not receive jet %d", i));
1333       continue;
1334     }
1335
1336     if (!AcceptJet(jet1))
1337       continue;
1338
1339     if (!AcceptBiasJet(jet1))
1340       continue;
1341
1342     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
1343       continue;
1344
1345     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
1346       continue;
1347
1348     if (!jet1->MatchedJet()) {
1349       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1350       if (!fRhoName.IsNull())
1351         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1352     }
1353
1354     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1355     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
1356
1357     if (naccJets1 < fNLeadingJets)
1358       fHistLeadingJets1PtArea->Fill(jet1->Area(), jet1->Pt());
1359
1360     if (!fRhoName.IsNull()) {
1361       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1362
1363       if (naccJets1 < fNLeadingJets)
1364         fHistLeadingJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
1365     }
1366
1367     if (fTracks) {
1368       for (Int_t it = 0; it < jet1->GetNumberOfTracks(); it++) {
1369         AliVParticle *track1 = jet1->TrackAt(it, fTracks2);
1370         if (track1)
1371           fHistJets1ZvsPt->Fill(track1->Pt() / jet1->Pt(), jet1->Pt());
1372       }
1373     }
1374
1375     if (fCaloClusters) {
1376       for (Int_t ic = 0; ic < jet1->GetNumberOfClusters(); ic++) {
1377         AliVCluster *cluster1 = jet1->ClusterAt(ic, fCaloClusters);
1378
1379         if (cluster1) {
1380           TLorentzVector nPart1;
1381           cluster1->GetMomentum(nPart1, fVertex);
1382           fHistJets2ZvsPt->Fill(nPart1.Et() / jet1->Pt(), jet1->Pt());
1383         }
1384       }
1385     }
1386
1387     fHistJets1NEFvsPt->Fill(jet1->NEF(), jet1->Pt());
1388     fHistJets1CEFvsCEFPt->Fill(1-jet1->NEF(), (1-jet1->NEF())*jet1->Pt());
1389
1390     naccJets1++;
1391   }
1392
1393   return kTRUE;
1394 }