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