PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx

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