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