PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetEmbed.cxx

   1 #include <TCanvas.h>
   2 #include <TChain.h>
   3 #include <TFormula.h>
   4 #include <TH1.h>
   5 #include <TH2.h>
   6 #include <TH3.h>
   7 #include <TProfile2D.h>
   8 #include <THnSparse.h>
   9 #include <TROOT.h>
  10 #include <TTree.h>
  11 #include <TArrayI.h>
  12 #include <TClonesArray.h>
  13 #include <TRandom3.h>
  14 #include <TFile.h>
  15 #include <TF1.h>
  16 #include <TLorentzVector.h>
  17 #include <TParameter.h>
  18
  19 #include "AliAODEvent.h"
  20 #include "AliAODInputHandler.h"
  21 #include "AliAnalysisManager.h"
  22 #include "AliAnalysisTask.h"
  23 #include "AliCentrality.h"
  24 #include "AliAnalysisTaskHJetEmbed.h"
  25 #include "AliESDEvent.h"
  26 #include "AliESDInputHandler.h"
  27 #include "AliVParticle.h"
  28 #include "AliVTrack.h"
  29 #include "AliInputEventHandler.h"
  30 #include "AliMCEvent.h"
  31 #include "AliStack.h"
  32 #include "AliGenEventHeader.h"
  33 #include "AliGenPythiaEventHeader.h"
  34 #include "AliLog.h"
  35 #include "AliRhoParameter.h"
  36 #include "AliNamedString.h"
  37
  38 #include "AliEmcalJet.h"
  39
  40 #include <iostream>
  41 using std::cout;
  42 using std::cerr;
  43 using std::endl;
  44
  45 ClassImp(AliAnalysisTaskHJetEmbed)
  46
  47 const Double_t pi = TMath::Pi();
  48 //const Double_t areaCut[4] = {0.1, 0.23, 0.4, 0.63};
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed() :
  52   AliAnalysisTaskSE(),
  53   fVerbosity(0), fAnaType(1), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
  54   fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
  55   fMCParticleArrName(""), fMCParticleArray(0x0),
  56   fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1),
  57   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
  58   fTTtype(0),
  59   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
  60   fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
  61   fRhoName(""), fRho(0x0), fRhoValue(0),
  62   fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
  63   fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
  64   fRunPL(kFALSE), fRunDL(kTRUE),
  65   fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
  66 {
  67   // Constructor
  68
  69   // Define input and output slots here
  70   // Input slot #0 works with a TChain
  71   //DefineInput(0, TChain::Class());
  72   //DefineOutput(1, TList::Class());
  73   // Output slot #0 id reserved by the base class for AOD
  74
  75   for(Int_t i=0; i<kNTrig; i++)
  76     {
  77       fhVtxZ[i]                  = 0x0;
  78       fhCentrality[i]            = 0x0;
  79       fhRhoVsCent[i]             = 0x0;
  80
  81       fhDLJetPtVsCent[i]         = 0x0;
  82       fhPLJetPtVsCent[i]         = 0x0;
  83
  84       fhPLTT[i]                  = 0x0;
  85       fhDLTT[i]                   = 0x0;
  86       fhPLHJet[i]                = 0x0;
  87       fhDLHJet[i]                = 0x0;
  88       fhTTPtQA[i]                = 0x0;
  89       fhTTPt[i]                  = 0x0;
  90       fhHJet[i]                  = 0x0;
  91
  92       fhJetPtGeoMatch[i]         = 0x0;
  93       fhJetPtEnMatch[i]          = 0x0;
  94       fhJetPhiGeoMatch[i]        = 0x0;
  95       fhJetPhiEnMatch[i]         = 0x0;
  96     }
  97   for(Int_t i=0; i<kNTT; i++)
  98     {
  99       if(i==0)
 100         { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
 101       else
 102         { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
 103     }
 104 }
 105
 106 //________________________________________________________________________
 107 AliAnalysisTaskHJetEmbed::AliAnalysisTaskHJetEmbed(const char *name) :
 108   AliAnalysisTaskSE(name),
 109   fVerbosity(0), fAnaType(1), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
 110   fEvent(0), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
 111   fMCParticleArrName(""), fMCParticleArray(0x0),
 112   fTrackArrName(""), fTrackArray(0x0), fTriggerTrkIndex(-1),
 113   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
 114   fTTtype(0),
 115   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""),
 116   fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
 117   fRhoName(""), fRho(0x0), fRhoValue(0),
 118   fPtHardBinName(0x0), fPtHardBin(-1), fRandom(0x0),
 119   fRunQA(kTRUE), fRunHJet(kTRUE), fRunMatch(kTRUE),
 120   fRunPL(kFALSE), fRunDL(kTRUE),
 121   fOutputList(0x0), fhEventStat(0x0), fhPtHardBins(0x0)
 122 {
 123   // Constructor
 124
 125   DefineInput(0, TChain::Class());
 126   DefineOutput(1, TList::Class());
 127
 128   for(Int_t i=0; i<kNTrig; i++)
 129     {
 130       fhVtxZ[i]                  = 0x0;
 131       fhCentrality[i]            = 0x0;
 132       fhRhoVsCent[i]             = 0x0;
 133
 134       fhDLJetPtVsCent[i]         = 0x0;
 135       fhPLJetPtVsCent[i]         = 0x0;
 136
 137       fhPLTT[i]                  = 0x0;
 138       fhDLTT[i]                   = 0x0;
 139       fhPLHJet[i]                = 0x0;
 140       fhDLHJet[i]                = 0x0;
 141       fhTTPtQA[i]                = 0x0;
 142       fhTTPt[i]                  = 0x0;
 143       fhHJet[i]                  = 0x0;
 144
 145       fhJetPtGeoMatch[i]         = 0x0;
 146       fhJetPtEnMatch[i]          = 0x0;
 147       fhJetPhiGeoMatch[i]        = 0x0;
 148       fhJetPhiEnMatch[i]         = 0x0;
 149     }
 150   for(Int_t i=0; i<kNTT; i++)
 151     {
 152       if(i==0)
 153         { fMinTTPt[i] = 19; fMaxTTPt[i] = 25; }
 154       else
 155         { fMinTTPt[i] = -1; fMaxTTPt[i] = -1; }
 156     }
 157 }
 158 //________________________________________________________________________
 159 AliAnalysisTaskHJetEmbed::~AliAnalysisTaskHJetEmbed()
 160 {
 161   //Destructor
 162   if(fRho)         delete fRho;
 163   if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
 164 }
 165
 166 //________________________________________________________________________
 167 void AliAnalysisTaskHJetEmbed::UserCreateOutputObjects()
 168 {
 169   // Create histograms
 170
 171   const Int_t nJetPtBins = 40;
 172   const Float_t lowJetPtBin=-50, upJetPtBin=150;
 173
 174   const Int_t nTrkPtBins = 100;
 175   const Float_t lowTrkPtBin=0, upTrkPtBin=100;
 176
 177
 178   // QA
 179   const Int_t dimJetqa = 3;
 180   const Int_t nBinsJetqa[dimJetqa]     = {nJetPtBins,  30, 11};
 181   const Double_t lowBinJetqa[dimJetqa] = {lowJetPtBin, 0,  0};
 182   const Double_t hiBinJetqa[dimJetqa]  = {upJetPtBin,  30, 11};
 183
 184   // h+jet
 185   const Int_t dimTT = 3;
 186   const Int_t nBinsTT[dimTT]     = {nTrkPtBins,  30, 11};
 187   const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0,  0};
 188   const Double_t hiBinTT[dimTT]  = {upTrkPtBin,  30, 11};
 189
 190   const Int_t dimHJet = 6;
 191   const Int_t nBinsHJet[dimHJet]     = {nTrkPtBins,  nJetPtBins,  140,     8,   30, 11};
 192   const Double_t lowBinHJet[dimHJet] = {lowTrkPtBin, lowJetPtBin, pi-4.95, 0,   0,  0};
 193   const Double_t hiBinHJet[dimHJet]  = {upTrkPtBin,  upJetPtBin,  pi+2.05, 0.8, 30, 11};
 194
 195   // Match
 196   const Int_t dimMthPt = 6;
 197   const Int_t nBinsMthPt[dimMthPt]     = {nJetPtBins,  nJetPtBins,  20,    20, 30, 11};
 198   const Double_t lowBinMthPt[dimMthPt] = {lowJetPtBin, lowJetPtBin, -0.95, 0,  0,  0};
 199   const Double_t hiBinMthPt[dimMthPt]  = {upJetPtBin,  upJetPtBin,  1.05,  1,  30, 11};
 200
 201   const Int_t dimMthPhi = 7;
 202   const Int_t nBinsMthPhi[dimMthPhi]     = {nTrkPtBins,  nJetPtBins/2,  70,      70,      10, 1,  11};
 203   const Double_t lowBinMthPhi[dimMthPhi] = {lowTrkPtBin, lowJetPtBin,   pi-4.95, pi-4.95, 0,  0,  0};
 204   const Double_t hiBinMthPhi[dimMthPhi]  = {upTrkPtBin,  upJetPtBin,    pi+2.05, pi+2.05, 1,  10, 11};
 205
 206   OpenFile(1);
 207   fOutputList = new TList();
 208   fOutputList->SetOwner(1);
 209
 210   fhEventStat = new TH1F("fhEventStat","Event statistics for jet analysis",9,0,9);
 211   fhEventStat->GetXaxis()->SetBinLabel(1,"All");
 212   fhEventStat->GetXaxis()->SetBinLabel(2,"PS");
 213   fhEventStat->GetXaxis()->SetBinLabel(3,"Vtx");
 214   fhEventStat->GetXaxis()->SetBinLabel(4,"Vtx+10cm");
 215   fhEventStat->GetXaxis()->SetBinLabel(5,"kMB");
 216   fhEventStat->GetXaxis()->SetBinLabel(6,"kCentral");
 217   fhEventStat->GetXaxis()->SetBinLabel(7,"kSemiCentral");
 218   fhEventStat->GetXaxis()->SetBinLabel(8,"kEMCEGA");
 219   fhEventStat->GetXaxis()->SetBinLabel(9,"kEMCEJE");
 220   fOutputList->Add(fhEventStat);
 221
 222   fhPtHardBins = new TH1F("fhPtHardBins","Number of events in each pT hard bin",11,0,11);
 223   fOutputList->Add(fhPtHardBins);
 224
 225   const char *triggerName[kNTrig] = {"kMB","kEGA","kEJE"};
 226
 227   for(Int_t i=0; i<kNTrig; i++)
 228     {
 229       if( i!=0 ) continue;
 230
 231       fhVtxZ[i] = new TH1F(Form("%s_fhVtxZ",triggerName[i]),Form("%s: z distribution of event vertexz;z(cm)",triggerName[i]),400,-20,20);
 232       fOutputList->Add(fhVtxZ[i]);
 233
 234       fhCentrality[i] = new TH1F(Form("%s_fhCentrality",triggerName[i]),Form("%s: Event centrality;centrality",triggerName[i]),100,0,100);
 235       fOutputList->Add(fhCentrality[i]);
 236
 237       fhRhoVsCent[i] = new TH2F(Form("%s_fhRhoVsCent",triggerName[i]),Form("%s: Rho vs centrality (R=%1.1f);centrality;Rho",triggerName[i],fRadius),100,0,100,300,0,300);
 238       fOutputList->Add(fhRhoVsCent[i]);
 239
 240       if(fRunQA)
 241         {
 242           fhPLJetPtVsCent[i] = new THnSparseF(Form("%s_fhPLJetPtVsCent",triggerName[i]),Form("PYTHIA: jet p_{T} vs centrality vs pT hard bin (particle-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
 243           fOutputList->Add(fhPLJetPtVsCent[i]);
 244
 245           fhDLJetPtVsCent[i] = new THnSparseF(Form("%s_fhDLJetPtVsCent",triggerName[i]),Form("PYTHIA: jet p_{T} vs centrality vs pT hard bin (detector-level,R=%1.1f);p_{T,jet}^{ch} (GeV/c);centrality;pT hard bin",fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
 246           fOutputList->Add(fhDLJetPtVsCent[i]);
 247         }
 248
 249       if(fRunHJet)
 250         {
 251           fhPLTT[i] = new THnSparseF(Form("%s_fhPLTT",triggerName[i]),Form("PYTHIA: TT p_{T} vs centrality vs pT hard bin (particle-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
 252           fOutputList->Add(fhPLTT[i]);
 253
 254           fhDLTT[i] = new THnSparseF(Form("%s_fhDLTT",triggerName[i]),Form("PYTHIA: TT p_{T} vs centrality vs pT hard bin (detector-level);p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
 255           fOutputList->Add(fhDLTT[i]);
 256
 257           fhTTPt[i] = new THnSparseF(Form("%s_fhTTPt",triggerName[i]),Form("Embedded: TT p_{T} vs centrality vs pT hard bin;p_{T,TT}^{ch} (GeV/c);centrality;pT hard bin"),dimTT,nBinsTT,lowBinTT,hiBinTT);
 258           fOutputList->Add(fhTTPt[i]);
 259
 260           fhPLHJet[i] = new THnSparseF(Form("%s_fhPLHJet",triggerName[i]),Form("PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (particle-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
 261           fOutputList->Add(fhPLHJet[i]);
 262
 263           fhDLHJet[i] = new THnSparseF(Form("%s_fhDLHJet",triggerName[i]),Form("PYTHIA: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (detector-level, R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
 264           fOutputList->Add(fhDLHJet[i]);
 265
 266           fhHJet[i] = new THnSparseF(Form("%s_fhHJet",triggerName[i]),Form("Embedded: TT p_{T} vs jet p_{T} vs #Delta#varphi vs jet area vs centrality vs pT hard bin (R=%1.1f);p_{T,TT}^{ch} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;Area;centrality;pT hard bin",fRadius),dimHJet,nBinsHJet,lowBinHJet,hiBinHJet);
 267           fOutputList->Add(fhHJet[i]);
 268         }
 269
 270       if(fRunMatch)
 271         {
 272           fhJetPtGeoMatch[i] = new THnSparseF(Form("%s_fhJetPtGeoMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
 273           fOutputList->Add(fhJetPtGeoMatch[i]);
 274
 275           fhJetPtEnMatch[i] = new THnSparseF(Form("%s_fhJetPtEnMatch",triggerName[i]),Form("Embed: generated p_{T,jet} vs reconstructed p_{T,jet} vs jet p_{T} difference vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c);(p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{gen};dR;centrality;pT hard bin",fRadius),dimMthPt,nBinsMthPt,lowBinMthPt,hiBinMthPt);
 276           fOutputList->Add(fhJetPtEnMatch[i]);
 277
 278           fhJetPhiGeoMatch[i] = new THnSparseF(Form("%s_fhJetPhiGeoMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
 279           fOutputList->Add(fhJetPhiGeoMatch[i]);
 280
 281           fhJetPhiEnMatch[i] = new THnSparseF(Form("%s_fhJetPhiEnMatch",triggerName[i]),Form("Embed: p_{T,TT} vs p_{T,jet}^{det} vs #Delta#varphi_{TT} vs #Delta#varphi_{jet} vs dR vs centrality vs pT hard bin (R=%1.1f);p_{T,TT} (GeV/c);p_{T,jet}^{det} (GeV/c);#Delta#varphi_{TT};#Delta#varphi_{jet};dR;centrality;pThard bin",fRadius),dimMthPhi,nBinsMthPhi,lowBinMthPhi,hiBinMthPhi);
 282           fOutputList->Add(fhJetPhiEnMatch[i]);
 283         }
 284     }
 285
 286   //error calculation in THnSparse
 287   Int_t nObj = fOutputList->GetEntries();
 288   for(Int_t i=0; i<nObj; i++)
 289     {
 290       TObject *obj = (TObject*) fOutputList->At(i);
 291       if (obj->IsA()->InheritsFrom( "THnSparse" ))
 292         {
 293           THnSparseF *hn = (THnSparseF*)obj;
 294           hn->Sumw2();
 295         }
 296     }
 297
 298   fRandom = new TRandom3();
 299
 300   PrintConfig();
 301   PostData(1, fOutputList);
 302 }
 303
 304 //________________________________________________________________________
 305 void AliAnalysisTaskHJetEmbed::UserExec(Option_t *)
 306 {
 307   // Main loop, called for each event.
 308
 309   AliDebug(5,"Entering UserExec");
 310   fTriggerType = -1;
 311   fEvent = InputEvent();
 312   if (!fEvent)
 313     {
 314       AliError("Input event not available");
 315       return;
 316     }
 317   AliDebug(5,"Got the input event");
 318   if(!fPtHardBinName)
 319     {
 320       // Get embedded pt hard bin number
 321       fPtHardBinName = static_cast<AliNamedString*>(fEvent->FindListObject("AODEmbeddingFile"));
 322       if(!fPtHardBinName)
 323         {
 324           AliError("The object for pt hard bin information is not available!");
 325           return;
 326         }
 327     }
 328   TString fileName = fPtHardBinName->GetString();
 329   fileName.Remove(0,50);
 330   fileName.Remove(fileName.Index("/"));
 331   fPtHardBin = fileName.Atoi();
 332   fhEventStat->Fill(0.5);
 333   UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 334
 335   if(fPeriod.Contains("lhc11h",TString::kIgnoreCase))
 336     {
 337       if (trigger & AliVEvent::kAnyINT)      { fTriggerType=0; }
 338       else if (trigger & AliVEvent::kCentral)     { fTriggerType=0; }
 339       else if (trigger & AliVEvent::kSemiCentral) { fTriggerType=0; }
 340       else if (trigger & AliVEvent::kEMCEGA)      { fTriggerType=1; }
 341       else if (trigger & AliVEvent::kEMCEJE)      { fTriggerType=2; }
 342     }
 343   else if(fPeriod.Contains("lhc10h",TString::kIgnoreCase))
 344     {
 345       if (trigger & AliVEvent::kAnyINT)   { fTriggerType=0; }
 346     }
 347   else if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase))
 348     {
 349       fTriggerType=0;
 350     }
 351
 352   if(fTriggerType==-1) return;
 353   fhEventStat->Fill(1.5);
 354
 355   // Vertex cut
 356   const AliVVertex* vtx = fEvent->GetPrimaryVertex();
 357   if (!vtx || vtx->GetNContributors()<1) return;
 358   fhEventStat->Fill(2.5);
 359   fhVtxZ[fTriggerType]->Fill(vtx->GetZ());
 360   if (TMath::Abs(vtx->GetZ())>fMaxVtxZ)  return;
 361   fhEventStat->Fill(3.5);
 362
 363   if(fTriggerType==0)
 364     {
 365       if(trigger & AliVEvent::kCentral) fhEventStat->Fill(5.5);
 366       else if (trigger & AliVEvent::kCentral) fhEventStat->Fill(6.5);
 367       else fhEventStat->Fill(4.5);
 368     }
 369   if(fTriggerType==1) fhEventStat->Fill(7.5);
 370   if(fTriggerType==2) fhEventStat->Fill(8.5);
 371
 372   // GetCentrality
 373   if(fCollisionSystem=="PbPb")
 374     {
 375       AliCentrality *centrality = fEvent->GetCentrality();
 376       if (centrality)
 377         fCentrality = centrality->GetCentralityPercentile("V0M");
 378       else
 379         fCentrality = 99;
 380     }
 381   else if(fCollisionSystem=="pp")
 382     {
 383       fCentrality = 0;
 384     }
 385
 386   // Get track collection and run QA
 387   if (!fTrackArray)
 388     {
 389       fTrackArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTrackArrName));
 390       if (!fTrackArray)
 391         {
 392           AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
 393           return;
 394         }
 395       if (!fTrackArray->GetClass()->GetBaseClass("AliVParticle"))
 396         {
 397           AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrName.Data()));
 398           fTrackArray = 0;
 399           return;
 400         }
 401     }
 402
 403   // Get MC particle array
 404   if (fRunPL && !fMCParticleArray)
 405     {
 406       fMCParticleArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fMCParticleArrName));
 407       if (!fMCParticleArray)
 408         {
 409           AliError(Form("Could not retrieve tracks %s!", fMCParticleArrName.Data()));
 410           return;
 411         }
 412     }
 413
 414   // Get Rho value
 415   if(fCollisionSystem=="PbPb")
 416     {
 417       if(!fRho && !fRhoName.IsNull())
 418         {
 419           fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
 420           if(!fRho)
 421             {
 422               AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
 423               return;
 424             }
 425         }
 426       if(fRho) fRhoValue = fRho->GetVal();
 427     }
 428   else if(fCollisionSystem=="pp")
 429     {
 430       fRhoValue = 0;
 431     }
 432
 433   // Get jet collection
 434   if (!fJetArray && !fJetArrName.IsNull())
 435     {
 436       fJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fJetArrName));
 437       if (!fJetArray)
 438         {
 439           AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
 440           return;
 441         }
 442       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 443         {
 444           AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrName.Data()));
 445           fJetArray = 0;
 446           return;
 447         }
 448     }
 449   // Get particle-level jet array
 450   if (fRunPL && !fPLJetArray && !fPLJetArrName.IsNull())
 451     {
 452       fPLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fPLJetArrName));
 453       if (!fPLJetArray)
 454         {
 455           AliError(Form("%s: Could not retrieve jets %s!", GetName(), fPLJetArrName.Data()));
 456           return;
 457         }
 458       if (!fPLJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 459         {
 460           AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fPLJetArrName.Data()));
 461           fPLJetArray = 0;
 462           return;
 463         }
 464     }
 465  // Get detector-level jet array
 466   if (fRunDL && !fDLJetArray && !fDLJetArrName.IsNull())
 467     {
 468       fDLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fDLJetArrName));
 469       if (!fDLJetArray)
 470         {
 471           AliError(Form("%s: Could not retrieve jets %s!", GetName(), fDLJetArrName.Data()));
 472           return;
 473         }
 474       if (!fDLJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 475         {
 476           AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fDLJetArrName.Data()));
 477           fDLJetArray = 0;
 478           return;
 479         }
 480     }
 481
 482   fhCentrality[fTriggerType]->Fill(fCentrality);
 483   fhRhoVsCent[fTriggerType]->Fill(fCentrality,fRhoValue);
 484   fhPtHardBins->Fill(fPtHardBin);
 485
 486   if(fRunQA) RunQA();
 487   if(fRunHJet)
 488     {
 489       for(Int_t i=0; i<kNTT; i++)
 490         RunHJet(fMinTTPt[i],fMaxTTPt[i]);
 491     }
 492
 493   PostData(1, fOutputList);
 494   return;
 495 }
 496
 497
 498 //________________________________________________________________________
 499 void AliAnalysisTaskHJetEmbed::RunHJet(const Double_t minPt, const Double_t maxPt)
 500 {
 501   TArrayI arr;
 502   Int_t counter = 0;
 503   Int_t indexPL = -1;
 504
 505   if(fRunPL)
 506     {
 507       // Find trigger track on particle level
 508       const Int_t nParticles = fMCParticleArray->GetEntries();
 509       Double_t maxPLPt = -1;
 510       arr.Set(nParticles);
 511       counter = 0;
 512       for(Int_t iPart=0; iPart<nParticles; iPart++)
 513         {
 514           AliVParticle *t = static_cast<AliVParticle*>(fMCParticleArray->At(iPart));
 515           //if(!t || t->Charge()==0) continue;
 516           //if(!AcceptTrack(t)) continue;
 517           Double_t pt = t->Pt();
 518           if(fTTtype==0) // single inclusive triggers
 519             {
 520               if (pt<maxPt && pt>=minPt)
 521                 {
 522                   arr.AddAt(iPart,counter);
 523                   counter++;
 524                 }
 525             }
 526           else if(fTTtype==1) // leading triggers
 527             {
 528               if(maxPLPt<pt)
 529                 {
 530                   maxPLPt = pt;
 531                   indexPL = iPart;
 532                 }
 533             }
 534         }
 535       arr.Set(counter);
 536       if(fTTtype==0)
 537         {
 538           if(counter==0) indexPL = -1;
 539           else if(counter==1) indexPL = arr.At(0);
 540           else
 541             {
 542               Double_t pro = fRandom->Uniform() * counter;
 543               indexPL = arr.At(TMath::FloorNint(pro));
 544             }
 545         }
 546       arr.Reset();
 547     }
 548
 549
 550   // Find trigger track on detector level and after embedding
 551   const Int_t Ntracks = fTrackArray->GetEntries();
 552   Double_t maxDLPt = 0;
 553   Int_t indexDL = -1;
 554   arr.Set(Ntracks);
 555   counter = 0;
 556   for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks)
 557     {
 558       AliVParticle *t = static_cast<AliVParticle*>(fTrackArray->At(iTracks));
 559       if(!t || t->Charge()==0) continue;
 560       if(!AcceptTrack(t)) continue;
 561       if(t->GetLabel()!=0)
 562         {
 563           //cout<<iTracks<<"  "<<t->Pt()<<" "<<t->GetLabel()<<endl;
 564           Double_t pt = t->Pt();
 565           if(fTTtype==0)
 566             {
 567               if (pt<maxPt && pt>=minPt)
 568                 {
 569                   arr.AddAt(iTracks,counter);
 570                   counter++;
 571                 }
 572             }
 573           else if(fTTtype==1)
 574             {
 575               if(maxDLPt<pt)
 576                 {
 577                   maxDLPt = pt;
 578                   indexDL = iTracks;
 579                 }
 580             }
 581         }
 582     }
 583   arr.Set(counter);
 584   if(fTTtype==0)
 585     {
 586       if(counter==0) indexDL = -1;
 587       else if(counter==1) indexDL = arr.At(0);
 588       else
 589         {
 590           Double_t pro = fRandom->Uniform() * counter;
 591           indexDL = arr.At(TMath::FloorNint(pro));
 592         }
 593     }
 594   arr.Reset();
 595
 596   AliDebug(2,Form("TT indices: PL=%d, DL=%d\n",indexPL,indexDL));
 597
 598   // Run h+jet
 599   if(fRunPL)  FillHJetCor(fMCParticleArray, indexPL, fPLJetArray, fhPLTT[fTriggerType], fhPLHJet[fTriggerType], kFALSE);
 600   if(fRunDL)  FillHJetCor(fTrackArray,      indexDL, fDLJetArray, fhDLTT[fTriggerType], fhDLHJet[fTriggerType], kFALSE);
 601   FillHJetCor(fTrackArray,      indexDL, fJetArray,   fhTTPt[fTriggerType], fhHJet[fTriggerType],   kTRUE);
 602
 603   if(fRunMatch) RunMatch(fTrackArray, indexDL);
 604 }
 605
 606 //________________________________________________________________________
 607 void  AliAnalysisTaskHJetEmbed::FillHJetCor(const TClonesArray *tracks, const Int_t leadingIndex, const TClonesArray *jetArray, THnSparse *hTT, THnSparse *hn, Bool_t isBkg)
 608 {
 609   if(leadingIndex<0) return;
 610
 611   AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
 612   Double_t triggerPt = tt->Pt();
 613   Double_t fill1[] = {triggerPt, fCentrality, static_cast<Double_t>(fPtHardBin) };
 614   hTT->Fill(fill1);
 615   AliDebug(2,Form("Found a trigger with pt = %2.2f",triggerPt));
 616
 617   Double_t triggerPhi = tt->Phi();
 618   if(triggerPhi<0) triggerPhi += 2*pi;
 619   Int_t nJets = jetArray->GetEntries();
 620   for(Int_t ij=0; ij<nJets; ij++)
 621     {
 622       AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
 623       if(!jet) continue;
 624       if(!IsGoodJet(jet)) continue; // eta cut
 625       Double_t jetPhi = jet->Phi();
 626       Double_t jetPt  = jet->Pt();
 627       Double_t jetArea = jet->Area();
 628       Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
 629       Double_t fill[] = {triggerPt,jetPt-jetArea*fRhoValue,dPhi,jetArea,fCentrality,static_cast<Double_t>(fPtHardBin)};
 630       if(!isBkg) fill[1] = jetPt;
 631       AliDebug(10,"Fill the histograms");
 632       hn->Fill(fill);
 633     }
 634 }
 635
 636 //________________________________________________________________________
 637 void AliAnalysisTaskHJetEmbed::RunMatch(const TClonesArray *tracks, const Int_t leadingIndex)
 638 {
 639   if(leadingIndex<0) return;
 640
 641   if(!fDLJetArray || !fJetArray)
 642     {
 643       AliWarning("Jet array is not available.");
 644       return;
 645     }
 646
 647   AliVParticle *tt = (AliVParticle*) tracks->At(leadingIndex);
 648   Double_t dR = 999, fraction = -1;
 649   Int_t nJets = fDLJetArray->GetEntries();
 650   for(Int_t ij=0; ij<nJets; ij++)
 651     {
 652       AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
 653       if(!jet) continue;
 654       if(!IsGoodJet(jet)) continue; // eta cut
 655       Double_t jetPt = jet->Pt();
 656       if(jetPt<10) continue;
 657
 658       // energy matching
 659       Int_t mthJetIndexEn = FindEnergyMatchedJet(jet,fJetArray,dR,fraction);
 660       if(mthJetIndexEn>-1 && fraction>0.5)
 661         {
 662           AliEmcalJet* jetMthEn = dynamic_cast<AliEmcalJet*>(fJetArray->At(mthJetIndexEn));
 663           if(jetMthEn)
 664             {
 665               Double_t fill[] = {tt->Pt(),jetPt,CalculateDPhi(tt->Phi(),jet->Phi()),CalculateDPhi(jetMthEn->Phi(),jet->Phi()),dR,fCentrality,static_cast<Double_t>(fPtHardBin)};
 666               fhJetPhiEnMatch[fTriggerType]->Fill(fill);
 667             }
 668         }
 669     }
 670
 671 }
 672
 673 //________________________________________________________________________
 674 Int_t AliAnalysisTaskHJetEmbed::FindGeoMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR)
 675 {
 676   dR = 999;
 677   if(!jetArray) return -1;
 678
 679   Int_t index = -1;
 680   Int_t nJets = jetArray->GetEntries();
 681   Double_t dRMax = 1;
 682   for(Int_t ij=0; ij<nJets; ij++)
 683     {
 684       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
 685       if(!jetTmp) continue;
 686       if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
 687       Double_t dPhi = GetDPhi(jet->Phi(),jetTmp->Phi());
 688       Double_t dEta = jet->Eta()-jetTmp->Eta();
 689       Double_t dRTmp = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 690       if(dRTmp<dRMax)
 691         {
 692           dRMax = dRTmp;
 693           index = ij;
 694         }
 695     }
 696   dR = dRMax;
 697   return index;
 698 }
 699
 700 //________________________________________________________________________
 701 Int_t AliAnalysisTaskHJetEmbed::FindEnergyMatchedJet(const AliEmcalJet* jet, const TClonesArray *jetArray, Double_t &dR, Double_t &fraction)
 702 {
 703   dR = 999;
 704   fraction=-1;
 705   if(!jetArray || !jet) return -1;
 706
 707   Int_t index = -1;
 708   Int_t nJets = jetArray->GetEntries();
 709   Double_t maxFrac = 0;
 710   Int_t nJetC = (Int_t)jet->GetNumberOfConstituents();
 711   Double_t jetPt = jet->Pt();
 712   for(Int_t ij=0; ij<nJets; ij++)
 713     {
 714       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
 715       if(!jetTmp) continue;
 716       if(TMath::Abs(jetTmp->Eta())>1) continue; // Generous eta cut
 717       if(GetJetDistance(jet,jetTmp)>1) continue;
 718
 719       Int_t nc = (Int_t)jetTmp->GetNumberOfConstituents();
 720       Double_t sumPt = 0;
 721       for(Int_t ic=0; ic<nc; ic++)
 722         {
 723           for(Int_t ijc=0; ijc<nJetC; ijc++)
 724             {
 725               if(jetTmp->TrackAt(ic)==jet->TrackAt(ijc))
 726                 {
 727                   AliVParticle *part = (AliVParticle*)jet->TrackAt(ijc,fTrackArray);
 728                   sumPt += part->Pt();
 729                 }
 730             }
 731         }
 732       Double_t frac = sumPt/jetPt;
 733       if(frac>maxFrac)
 734         {
 735           maxFrac = frac;
 736           index = ij;
 737         }
 738     }
 739   fraction = maxFrac;
 740
 741   if(index>0)
 742     {
 743       AliEmcalJet* jetTmp = dynamic_cast<AliEmcalJet*>(jetArray->At(index));
 744       if(jetTmp)
 745         dR = GetJetDistance(jet,jetTmp);
 746     }
 747   return index;
 748 }
 749
 750 //________________________________________________________________________
 751 Double_t AliAnalysisTaskHJetEmbed::CalculateDPhi(const Double_t phi1, const Double_t phi2)
 752 {
 753   Double_t dPhi = phi1-phi2;
 754   if(dPhi>2*pi)  dPhi -= 2*pi;
 755   if(dPhi<-2*pi) dPhi += 2*pi;
 756   if(dPhi<-0.5*pi) dPhi += 2*pi;
 757   if(dPhi>1.5*pi)  dPhi -= 2*pi;
 758   return dPhi;
 759 }
 760
 761 //________________________________________________________________________
 762 Double_t AliAnalysisTaskHJetEmbed::GetDPhi(const Double_t phi1, const Double_t phi2)
 763 {
 764   Double_t dPhi = TMath::Abs(phi1-phi2);
 765   if(dPhi>2*pi) dPhi -= 2*pi;
 766   if(dPhi>pi)   dPhi = 2*pi - dPhi;
 767   return dPhi;
 768 }
 769
 770 //________________________________________________________________________
 771 Double_t AliAnalysisTaskHJetEmbed::GetJetDistance(const AliEmcalJet *jet1, const AliEmcalJet* jet2)
 772 {
 773   Double_t dPhi = GetDPhi(jet1->Phi(),jet2->Phi());
 774   Double_t dEta = jet1->Eta()-jet2->Eta();
 775   return TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 776 }
 777
 778 //________________________________________________________________________
 779 void AliAnalysisTaskHJetEmbed::RunQA()
 780 {
 781   if(!fPLJetArray)
 782     {
 783       AliWarning(Form("Particle-level jet array is not available: %s\n",fPLJetArrName.Data()));
 784     }
 785   else
 786     {
 787       Int_t nPLJets = fPLJetArray->GetEntries();
 788       for(Int_t ij=0; ij<nPLJets; ij++)
 789         {
 790           AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fPLJetArray->At(ij));
 791           if(!jet) continue;
 792           if(!IsGoodJet(jet)) continue; // eta cut
 793           Double_t jetPt = jet->Pt();
 794           Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
 795           fhPLJetPtVsCent[fTriggerType]->Fill(fill);
 796           AliDebug(5, Form("PL jet %d has (pt,eta,phi) = (%2.2f,%2.2f,%2.2f)",ij,jetPt,jet->Eta(),jet->Phi()));
 797         }
 798     }
 799
 800   if(!fDLJetArray)
 801     {
 802       AliWarning(Form("Detector-level jet array is not available: %s\n",fDLJetArrName.Data()));
 803     }
 804   else
 805     {
 806       Int_t nDLJets = fDLJetArray->GetEntries();
 807       for(Int_t ij=0; ij<nDLJets; ij++)
 808         {
 809           AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
 810           if(!jet) continue;
 811           if(!IsGoodJet(jet)) continue; // eta cut
 812           Double_t jetPt = jet->Pt();
 813           Double_t fill[] = {jetPt, fCentrality, static_cast<Double_t>(fPtHardBin)};
 814           fhDLJetPtVsCent[fTriggerType]->Fill(fill);
 815           AliDebug(5, Form("DL jet %d has pt = %2.2f",ij,jetPt));
 816         }
 817     }
 818 }
 819
 820
 821 //________________________________________________________________________
 822 Bool_t AliAnalysisTaskHJetEmbed::AcceptTrack(const AliVParticle *track)
 823 {
 824   if(track->Pt()<fMinTrkPt || track->Pt()>fMaxTrkPt) return kFALSE;
 825   if(track->Eta()<fMinTrkEta || track->Eta()>fMaxTrkEta) return kFALSE;
 826   if(track->Phi()<fMinTrkPhi || track->Phi()>fMaxTrkPhi) return kFALSE;
 827   return kTRUE;
 828 }
 829
 830 //________________________________________________________________________
 831 Bool_t AliAnalysisTaskHJetEmbed::IsGoodJet(const AliEmcalJet* jet)
 832 {
 833   Double_t etaCut = (0.9-fRadius>0.5)?0.5:0.9-fRadius;
 834   if(TMath::Abs(jet->Eta())>etaCut) return kFALSE;
 835   return kTRUE;
 836 }
 837
 838 //
 839 //________________________________________________________________________
 840 //
 841 void AliAnalysisTaskHJetEmbed::PrintConfig()
 842 {
 843   const char *decision[2] = {"no","yes"};
 844   const char *TTtype[2] = {"Single inclusive","Leading"};
 845   printf("\n\n===== h-jet analysis configuration =====\n");
 846   printf("Input event type: %s - %s\n",fCollisionSystem.Data(),fPeriod.Data());
 847   printf("Track pt range: %2.2f < pt < %2.2f\n",fMinTrkPt, fMaxTrkPt);
 848   printf("Track eta range: %2.1f < eta < %2.1f\n",fMinTrkEta, fMaxTrkEta);
 849   printf("Track phi range: %2.0f < phi < %2.0f\n",fMinTrkPhi*TMath::RadToDeg(),fMaxTrkPhi*TMath::RadToDeg());
 850   printf("TT type: %s\n", TTtype[fTTtype]);
 851   for(Int_t i=0; i<kNTT; i++)
 852     printf("TT range %d:  %2.0f < pt < %2.0f\n", i+1, fMinTTPt[i], fMaxTTPt[i]);
 853   printf("Run QA: %s\n",decision[fRunQA]);
 854   printf("Run particle level: %s\n",decision[fRunPL]);
 855   printf("Run detector level: %s\n",decision[fRunDL]);
 856   printf("Run h+jet: %s\n",decision[fRunHJet]);
 857   printf("Run matching: %s\n",decision[fRunMatch]);
 858   printf("=======================================\n\n");
 859 }
 860
 861 //________________________________________________________________________
 862 Double_t AliAnalysisTaskHJetEmbed::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz)
 863 {
 864   return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
 865 }
 866
 867 //________________________________________________________________________
 868 void AliAnalysisTaskHJetEmbed::Terminate(Option_t *)
 869 {
 870   // Called once at the end of the query
 871 }