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