PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx

   1 // $Id$
   2 //
   3 // Emcal jet response matrix maker task.
   4 //
   5 // Author: S. Aiola
   6
   7 #include "AliJetResponseMaker.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TH3F.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliVTrack.h"
  17 #include "AliEmcalJet.h"
  18 #include "AliGenPythiaEventHeader.h"
  19 #include "AliLog.h"
  20 #include "AliMCEvent.h"
  21 #include "AliRhoParameter.h"
  22
  23 ClassImp(AliJetResponseMaker)
  24
  25 //________________________________________________________________________
  26 AliJetResponseMaker::AliJetResponseMaker() :
  27   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  28   fTracks2Name(""),
  29   fCalo2Name(""),
  30   fJets2Name(""),
  31   fRho2Name(""),
  32   fPtBiasJet2Track(0),
  33   fPtBiasJet2Clus(0),
  34   fAreCollections1MC(kFALSE),
  35   fAreCollections2MC(kTRUE),
  36   fMatching(kNoMatching),
  37   fMatchingPar(0),
  38   fJet2MinEta(-999),
  39   fJet2MaxEta(-999),
  40   fJet2MinPhi(-999),
  41   fJet2MaxPhi(-999),
  42   fSelectPtHardBin(-999),
  43   fPythiaHeader(0),
  44   fPtHardBin(0),
  45   fNTrials(0),
  46   fTracks2(0),
  47   fCaloClusters2(0),
  48   fJets2(0),
  49   fRho2(0),
  50   fRho2Val(0),
  51   fHistNTrials(0),
  52   fHistEvents(0),
  53   fHistJets1PhiEta(0),
  54   fHistJets1PtArea(0),
  55   fHistJets1CorrPtArea(0),
  56   fHistJets2PhiEta(0),
  57   fHistJets2PtArea(0),
  58   fHistJets2CorrPtArea(0),
  59   fHistMatchingLevelvsJet2Pt(0),
  60   fHistClosestDeltaEtaPhivsJet2Pt(0),
  61   fHistClosestDeltaPtvsJet2Pt(0),
  62   fHistClosestDeltaCorrPtvsJet2Pt(0),
  63   fHistNonMatchedJets1PtArea(0),
  64   fHistNonMatchedJets2PtArea(0),
  65   fHistNonMatchedJets1CorrPtArea(0),
  66   fHistNonMatchedJets2CorrPtArea(0),
  67   fHistJet1PtvsJet2Pt(0),
  68   fHistJet1CorrPtvsJet2CorrPt(0),
  69   fHistMissedJets2PtArea(0)
  70 {
  71   // Default constructor.
  72
  73   SetMakeGeneralHistograms(kTRUE);
  74 }
  75
  76 //________________________________________________________________________
  77 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
  78   AliAnalysisTaskEmcalJet(name, kTRUE),
  79   fTracks2Name("MCParticles"),
  80   fCalo2Name(""),
  81   fJets2Name("MCJets"),
  82   fRho2Name(""),
  83   fPtBiasJet2Track(0),
  84   fPtBiasJet2Clus(0),
  85   fAreCollections1MC(kFALSE),
  86   fAreCollections2MC(kTRUE),
  87   fMatching(kNoMatching),
  88   fMatchingPar(0.25),
  89   fJet2MinEta(-999),
  90   fJet2MaxEta(-999),
  91   fJet2MinPhi(-999),
  92   fJet2MaxPhi(-999),
  93   fSelectPtHardBin(-999),
  94   fPythiaHeader(0),
  95   fPtHardBin(0),
  96   fNTrials(0),
  97   fTracks2(0),
  98   fCaloClusters2(0),
  99   fJets2(0),
 100   fRho2(0),
 101   fRho2Val(0),
 102   fHistNTrials(0),
 103   fHistEvents(0),
 104   fHistJets1PhiEta(0),
 105   fHistJets1PtArea(0),
 106   fHistJets1CorrPtArea(0),
 107   fHistJets2PhiEta(0),
 108   fHistJets2PtArea(0),
 109   fHistJets2CorrPtArea(0),
 110   fHistMatchingLevelvsJet2Pt(0),
 111   fHistClosestDeltaEtaPhivsJet2Pt(0),
 112   fHistClosestDeltaPtvsJet2Pt(0),
 113   fHistClosestDeltaCorrPtvsJet2Pt(0),
 114   fHistNonMatchedJets1PtArea(0),
 115   fHistNonMatchedJets2PtArea(0),
 116   fHistNonMatchedJets1CorrPtArea(0),
 117   fHistNonMatchedJets2CorrPtArea(0),
 118   fHistJet1PtvsJet2Pt(0),
 119   fHistJet1CorrPtvsJet2CorrPt(0),
 120   fHistMissedJets2PtArea(0)
 121 {
 122   // Standard constructor.
 123
 124   SetMakeGeneralHistograms(kTRUE);
 125 }
 126
 127 //________________________________________________________________________
 128 AliJetResponseMaker::~AliJetResponseMaker()
 129 {
 130   // Destructor
 131 }
 132
 133 //________________________________________________________________________
 134 void AliJetResponseMaker::UserCreateOutputObjects()
 135 {
 136   // Create user objects.
 137
 138   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 139
 140   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 141   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 142
 143   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
 144   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 145   fHistNTrials->GetYaxis()->SetTitle("trials");
 146   fOutput->Add(fHistNTrials);
 147
 148   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 149   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 150   fHistEvents->GetYaxis()->SetTitle("total events");
 151   fOutput->Add(fHistEvents);
 152
 153   for (Int_t i = 1; i < 12; i++) {
 154     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 155     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 156   }
 157
 158   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 20, -2, 2, 32, 0, 6.4);
 159   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
 160   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
 161   fOutput->Add(fHistJets1PhiEta);
 162
 163   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 164   fHistJets1PtArea->GetXaxis()->SetTitle("area");
 165   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 166   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
 167   fOutput->Add(fHistJets1PtArea);
 168
 169   if (!fRhoName.IsNull()) {
 170     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 171     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
 172     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
 173     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 174     fOutput->Add(fHistJets1CorrPtArea);
 175   }
 176
 177   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 20, -2, 2, 32, 0, 6.4);
 178   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
 179   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
 180   fOutput->Add(fHistJets2PhiEta);
 181
 182   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 183   fHistJets2PtArea->GetXaxis()->SetTitle("area");
 184   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 185   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
 186   fOutput->Add(fHistJets2PtArea);
 187
 188   if (!fRho2Name.IsNull()) {
 189     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 190     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
 191     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 192     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 193     fOutput->Add(fHistJets2CorrPtArea);
 194   }
 195
 196   fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 197   fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
 198   fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 199   fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
 200   fOutput->Add(fHistMatchingLevelvsJet2Pt);
 201
 202   fHistClosestDeltaEtaPhivsJet2Pt = new TH3F("fHistClosestDeltaEtaPhivsJet2Pt", "fHistClosestDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
 203   fHistClosestDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
 204   fHistClosestDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
 205   fHistClosestDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
 206   fOutput->Add(fHistClosestDeltaEtaPhivsJet2Pt);
 207
 208   fHistClosestDeltaPtvsJet2Pt = new TH2F("fHistClosestDeltaPtvsJet2Pt", "fHistClosestDeltaPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 209   fHistClosestDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 210   fHistClosestDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 211   fHistClosestDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 212   fOutput->Add(fHistClosestDeltaPtvsJet2Pt);
 213
 214   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 215     fHistClosestDeltaCorrPtvsJet2Pt = new TH2F("fHistClosestDeltaCorrPtvsJet2Pt", "fHistClosestDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 216     fHistClosestDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 217     fHistClosestDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 218     fHistClosestDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 219     fOutput->Add(fHistClosestDeltaCorrPtvsJet2Pt);
 220   }
 221
 222   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 223   fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
 224   fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 225   fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
 226   fOutput->Add(fHistNonMatchedJets1PtArea);
 227
 228   fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 229   fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
 230   fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 231   fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
 232   fOutput->Add(fHistNonMatchedJets2PtArea);
 233
 234   if (!fRhoName.IsNull()) {
 235     fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 236     fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
 237     fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 238     fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 239     fOutput->Add(fHistNonMatchedJets1CorrPtArea);
 240   }
 241
 242   if (!fRho2Name.IsNull()) {
 243     fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 244     fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
 245     fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 246     fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 247     fOutput->Add(fHistNonMatchedJets2CorrPtArea);
 248   }
 249
 250   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 251   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
 252   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 253   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 254   fOutput->Add(fHistJet1PtvsJet2Pt);
 255
 256   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 257     if (fRhoName.IsNull())
 258       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 259     else if (fRho2Name.IsNull())
 260       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 261     else
 262       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 263     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 264     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
 265     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 266     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
 267   }
 268
 269   fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 270   fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
 271   fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 272   fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
 273   fOutput->Add(fHistMissedJets2PtArea);
 274
 275   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 276 }
 277
 278 //________________________________________________________________________
 279 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 280 {
 281   // Return true if jet is accepted.
 282
 283   if (jet->Pt() <= fJetPtCut)
 284     return kFALSE;
 285   if (jet->Area() <= fJetAreaCut)
 286     return kFALSE;
 287
 288   return kTRUE;
 289 }
 290
 291 //________________________________________________________________________
 292 Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
 293 {
 294   // Accept jet with a bias.
 295
 296   if (fLeadingHadronType == 0) {
 297     if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
 298   }
 299   else if (fLeadingHadronType == 1) {
 300     if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 301   }
 302   else {
 303     if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 304   }
 305
 306   return kTRUE;
 307 }
 308
 309 //________________________________________________________________________
 310 void AliJetResponseMaker::ExecOnce()
 311 {
 312   // Execute once.
 313
 314   if (!fJets2Name.IsNull() && !fJets2) {
 315     fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
 316     if (!fJets2) {
 317       AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
 318       return;
 319     }
 320     else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
 321       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
 322       fJets2 = 0;
 323       return;
 324     }
 325   }
 326
 327   if (!fTracks2Name.IsNull() && !fTracks2) {
 328     fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
 329     if (!fTracks2) {
 330       AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
 331       return;
 332     }
 333     else {
 334       TClass *cl = fTracks2->GetClass();
 335       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 336         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
 337         fTracks2 = 0;
 338         return;
 339       }
 340     }
 341   }
 342
 343   if (!fCalo2Name.IsNull() && !fCaloClusters2) {
 344     fCaloClusters2 =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 345     if (!fCaloClusters2) {
 346       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
 347       return;
 348     } else {
 349       TClass *cl = fCaloClusters2->GetClass();
 350       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
 351         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
 352         fCaloClusters2 = 0;
 353         return;
 354       }
 355     }
 356   }
 357
 358   if (!fRho2Name.IsNull() && !fRho2) {
 359     fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
 360     if (!fRho2) {
 361       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
 362       fInitialized = kFALSE;
 363       return;
 364     }
 365   }
 366
 367   if (fJet2MinEta == -999)
 368     fJet2MinEta = fJetMinEta - fJetRadius;
 369   if (fJet2MaxEta == -999)
 370     fJet2MaxEta = fJetMaxEta + fJetRadius;
 371   if (fJet2MinPhi == -999)
 372     fJet2MinPhi = fJetMinPhi - fJetRadius;
 373   if (fJet2MaxPhi == -999)
 374     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
 375
 376   AliAnalysisTaskEmcalJet::ExecOnce();
 377 }
 378
 379 //________________________________________________________________________
 380 Bool_t AliJetResponseMaker::IsEventSelected()
 381 {
 382   // Check if event is selected
 383
 384   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
 385     return kFALSE;
 386
 387   return AliAnalysisTaskEmcalJet::IsEventSelected();
 388 }
 389
 390 //________________________________________________________________________
 391 Bool_t AliJetResponseMaker::RetrieveEventObjects()
 392 {
 393   // Retrieve event objects.
 394
 395   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 396     return kFALSE;
 397
 398   if (fRho2)
 399     fRho2Val = fRho2->GetVal();
 400
 401   fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 402
 403   if (!fPythiaHeader)
 404     return kFALSE;
 405
 406   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 407   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 408
 409   Double_t pthard = fPythiaHeader->GetPtHard();
 410
 411   for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 412     if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
 413       break;
 414   }
 415
 416   fNTrials = fPythiaHeader->Trials();
 417
 418   return kTRUE;
 419 }
 420
 421 //________________________________________________________________________
 422 Bool_t AliJetResponseMaker::Run()
 423 {
 424   // Find the closest jets
 425
 426   if (fMatching == kNoMatching)
 427     return kTRUE;
 428
 429   DoJetLoop(fJets, fJets2, kFALSE);
 430   DoJetLoop(fJets2, fJets, kTRUE);
 431
 432   const Int_t nJets2 = fJets2->GetEntriesFast();
 433
 434   for (Int_t i = 0; i < nJets2; i++) {
 435
 436     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
 437
 438     if (!jet2) {
 439       AliError(Form("Could not receive jet %d", i));
 440       continue;
 441     }
 442
 443     if (!AcceptJet(jet2))
 444       continue;
 445
 446     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
 447       continue;
 448
 449     if (jet2->Pt() > fMaxBinPt)
 450       continue;
 451
 452     if (jet2->ClosestJet() && jet2->ClosestJet()->ClosestJet() == jet2 &&
 453         jet2->ClosestJetDistance() < fMatchingPar) {    // Matched jet found
 454       jet2->SetMatchedToClosest();
 455       jet2->ClosestJet()->SetMatchedToClosest();
 456     }
 457   }
 458
 459   return kTRUE;
 460 }
 461
 462 //________________________________________________________________________
 463 void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t first)
 464 {
 465   // Do the jet loop.
 466
 467   Int_t nJets1 = jets1->GetEntriesFast();
 468   Int_t nJets2 = jets2->GetEntriesFast();
 469
 470   for (Int_t i = 0; i < nJets1; i++) {
 471
 472     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 473
 474     if (!jet1) {
 475       AliError(Form("Could not receive jet %d", i));
 476       continue;
 477     }
 478
 479     if (!AcceptJet(jet1))
 480       continue;
 481
 482     if (first) {
 483      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 484         continue;
 485     }
 486     else {
 487       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 488         continue;
 489     }
 490
 491     for (Int_t j = 0; j < nJets2; j++) {
 492
 493       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
 494
 495       if (!jet2) {
 496         AliError(Form("Could not receive jet %d", j));
 497         continue;
 498       }
 499
 500       if (!AcceptJet(jet2))
 501         continue;
 502
 503       if (first) {
 504         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
 505           continue;
 506       }
 507       else {
 508         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 509           continue;
 510       }
 511
 512       Double_t d = GetMatchingLevel(jet1, jet2);
 513
 514       if (d < 0)
 515         continue;
 516
 517       if (d < jet1->ClosestJetDistance()) {
 518         jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
 519         jet1->SetClosestJet(jet2, d);
 520       }
 521       else if (d < jet1->SecondClosestJetDistance()) {
 522         jet1->SetSecondClosestJet(jet2, d);
 523       }
 524     }
 525   }
 526 }
 527
 528 //________________________________________________________________________
 529 Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2) const
 530 {
 531   Double_t r = -1;
 532
 533   switch (fMatching) {
 534   case kGeometrical:
 535     {
 536       Double_t deta = jet2->Eta() - jet1->Eta();
 537       Double_t dphi = jet2->Phi() - jet1->Phi();
 538       r = TMath::Sqrt(deta * deta + dphi * dphi);
 539     }
 540     break;
 541   case kMCLabel:
 542     AliError("MC label matching not implemented!");
 543     break;
 544   default:
 545     ;
 546   }
 547
 548   return r;
 549 }
 550
 551 //________________________________________________________________________
 552 Bool_t AliJetResponseMaker::FillHistograms()
 553 {
 554   // Fill histograms.
 555
 556   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
 557   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
 558
 559   const Int_t nJets2 = fJets2->GetEntriesFast();
 560
 561   for (Int_t i = 0; i < nJets2; i++) {
 562
 563     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
 564
 565     if (!jet2) {
 566       AliError(Form("Could not receive jet2 %d", i));
 567       continue;
 568     }
 569
 570     if (!AcceptJet(jet2))
 571       continue;
 572
 573     if (!AcceptBiasJet2(jet2))
 574       continue;
 575
 576     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
 577       continue;
 578
 579     if (jet2->Pt() > fMaxBinPt)
 580       continue;
 581
 582     if (jet2->MatchedJet()) {
 583
 584       if (!AcceptBiasJet(jet2->MatchedJet()) ||
 585           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
 586           jet2->MatchedJet()->Pt() > fMaxBinPt) {
 587         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 588       }
 589       else {
 590         fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
 591
 592         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
 593         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
 594         fHistClosestDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
 595
 596         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
 597         fHistClosestDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
 598
 599         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
 600
 601         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 602           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
 603           fHistClosestDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
 604           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
 605         }
 606       }
 607     }
 608     else {
 609       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 610       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 611
 612       if (!fRho2Name.IsNull())
 613         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
 614     }
 615
 616     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 617     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
 618
 619     if (!fRho2Name.IsNull())
 620       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
 621   }
 622
 623   const Int_t nJets1 = fJets->GetEntriesFast();
 624
 625   for (Int_t i = 0; i < nJets1; i++) {
 626
 627     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
 628
 629     if (!jet1) {
 630       AliError(Form("Could not receive jet %d", i));
 631       continue;
 632     }
 633
 634     if (!AcceptJet(jet1))
 635       continue;
 636
 637     if (!AcceptBiasJet(jet1))
 638       continue;
 639
 640     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
 641       continue;
 642
 643     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 644       continue;
 645
 646     if (!jet1->MatchedJet()) {
 647       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
 648       if (!fRhoName.IsNull())
 649         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
 650     }
 651
 652     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
 653     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
 654
 655     if (!fRhoName.IsNull())
 656       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
 657   }
 658
 659   return kTRUE;
 660 }