PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeDeriv.cxx

   1 //
   2 // Do subtraction for jet shapes using derivatives arXiv:1211:2811
   3 //
   4 // Author: M.Verweij
   5
   6 #include <TClonesArray.h>
   7 #include <TH1F.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THnSparse.h>
  11 #include <TF1.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TChain.h>
  16 #include <TSystem.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19 #include <TTree.h>
  20
  21 #include "AliVCluster.h"
  22 #include "AliVTrack.h"
  23 #include "AliEmcalJet.h"
  24 #include "AliRhoParameter.h"
  25 #include "AliLog.h"
  26 #include "AliEmcalParticle.h"
  27 #include "AliMCEvent.h"
  28 #include "AliAODEvent.h"
  29 #include "AliGenPythiaEventHeader.h"
  30 #include "AliAODMCHeader.h"
  31 #include "AliAnalysisManager.h"
  32 #include "AliJetContainer.h"
  33 #include "AliParticleContainer.h"
  34
  35 #include "AliAnalysisTaskJetShapeDeriv.h"
  36
  37 ClassImp(AliAnalysisTaskJetShapeDeriv)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv() :
  41   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeDeriv", kTRUE),
  42   fContainerBase(0),
  43   fContainerNoEmb(1),
  44   fMinFractionShared(0),
  45   fSingleTrackEmb(kFALSE),
  46   fCreateTree(kFALSE),
  47   fTreeJetBkg(),
  48   fJet1Vec(new TLorentzVector()),
  49   fJet2Vec(new TLorentzVector()),
  50   fArea(0),
  51   fAreaPhi(0),
  52   fAreaEta(0),
  53   fRho(0),
  54   fRhoM(0),
  55   fNConst(0),
  56   fM1st(0),
  57   fM2nd(0),
  58   fDeriv1st(0),
  59   fDeriv2nd(0),
  60   fMatch(0),
  61   fh2MSubMatch(0x0),
  62   fh2MSubPtRawAll(0x0),
  63   fh3MSubPtRawDRMatch(0x0),
  64   fh3MSubPtTrueDR(0x0),
  65   fh3MTruePtTrueDR(0x0),
  66   fh3PtTrueDeltaMDR(0x0),
  67   fh3PtTrueDeltaMRelDR(0x0),
  68   fhnMassResponse(0x0),
  69   fh2PtTrueSubFacV1(0x0),
  70   fh2PtRawSubFacV1(0x0),
  71   fh2PtCorrSubFacV1(0x0),
  72   fh2NConstSubFacV1(0x0),
  73   fh2PtTrueSubFacV2(0x0),
  74   fh2PtRawSubFacV2(0x0),
  75   fh2PtCorrSubFacV2(0x0),
  76   fh2NConstSubFacV2(0x0)
  77 {
  78   // Default constructor.
  79
  80   fh2MSubMatch         = new TH2F*[fNcentBins];
  81   fh2MSubPtRawAll      = new TH2F*[fNcentBins];
  82   fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
  83   fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
  84   fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
  85   fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
  86   fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
  87   fhnMassResponse      = new THnSparse*[fNcentBins];
  88   fh2PtTrueSubFacV1    = new TH2F*[fNcentBins];
  89   fh2PtRawSubFacV1     = new TH2F*[fNcentBins];
  90   fh2PtCorrSubFacV1    = new TH2F*[fNcentBins];
  91   fh2NConstSubFacV1    = new TH2F*[fNcentBins];
  92   fh2PtTrueSubFacV2    = new TH2F*[fNcentBins];
  93   fh2PtRawSubFacV2     = new TH2F*[fNcentBins];
  94   fh2PtCorrSubFacV2    = new TH2F*[fNcentBins];
  95   fh2NConstSubFacV2    = new TH2F*[fNcentBins];
  96
  97   for (Int_t i = 0; i < fNcentBins; i++) {
  98     fh2MSubMatch[i]         = 0;
  99     fh2MSubPtRawAll[i]      = 0;
 100     fh3MSubPtRawDRMatch[i]  = 0;
 101     fh3MSubPtTrueDR[i]      = 0;
 102     fh3MTruePtTrueDR[i]     = 0;
 103     fh3PtTrueDeltaMDR[i]    = 0;
 104     fh3PtTrueDeltaMRelDR[i] = 0;
 105     fhnMassResponse[i]      = 0;
 106     fh2PtTrueSubFacV1[i]    = 0;
 107     fh2PtRawSubFacV1[i]     = 0;
 108     fh2PtCorrSubFacV1[i]    = 0;
 109     fh2NConstSubFacV1[i]    = 0;
 110     fh2PtTrueSubFacV2[i]    = 0;
 111     fh2PtRawSubFacV2[i]     = 0;
 112     fh2PtCorrSubFacV2[i]    = 0;
 113     fh2NConstSubFacV2[i]    = 0;
 114   }
 115
 116   SetMakeGeneralHistograms(kTRUE);
 117   if(fCreateTree) DefineOutput(2, TTree::Class());
 118 }
 119
 120 //________________________________________________________________________
 121 AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv(const char *name) :
 122   AliAnalysisTaskEmcalJet(name, kTRUE),
 123   fContainerBase(0),
 124   fContainerNoEmb(1),
 125   fMinFractionShared(0),
 126   fSingleTrackEmb(kFALSE),
 127   fCreateTree(kFALSE),
 128   fTreeJetBkg(0),
 129   fJet1Vec(new TLorentzVector()),
 130   fJet2Vec(new TLorentzVector()),
 131   fArea(0),
 132   fAreaPhi(0),
 133   fAreaEta(0),
 134   fRho(0),
 135   fRhoM(0),
 136   fNConst(0),
 137   fM1st(0),
 138   fM2nd(0),
 139   fDeriv1st(0),
 140   fDeriv2nd(0),
 141   fMatch(0),
 142   fh2MSubMatch(0x0),
 143   fh2MSubPtRawAll(0x0),
 144   fh3MSubPtRawDRMatch(0x0),
 145   fh3MSubPtTrueDR(0x0),
 146   fh3MTruePtTrueDR(0x0),
 147   fh3PtTrueDeltaMDR(0x0),
 148   fh3PtTrueDeltaMRelDR(0x0),
 149   fhnMassResponse(0x0),
 150   fh2PtTrueSubFacV1(0x0),
 151   fh2PtRawSubFacV1(0x0),
 152   fh2PtCorrSubFacV1(0x0),
 153   fh2NConstSubFacV1(0x0),
 154   fh2PtTrueSubFacV2(0x0),
 155   fh2PtRawSubFacV2(0x0),
 156   fh2PtCorrSubFacV2(0x0),
 157   fh2NConstSubFacV2(0x0)
 158 {
 159   // Standard constructor.
 160
 161   fh2MSubMatch         = new TH2F*[fNcentBins];
 162   fh2MSubPtRawAll      = new TH2F*[fNcentBins];
 163   fh3MSubPtRawDRMatch  = new TH3F*[fNcentBins];
 164   fh3MSubPtTrueDR      = new TH3F*[fNcentBins];
 165   fh3MTruePtTrueDR     = new TH3F*[fNcentBins];
 166   fh3PtTrueDeltaMDR    = new TH3F*[fNcentBins];
 167   fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
 168   fhnMassResponse      = new THnSparse*[fNcentBins];
 169   fh2PtTrueSubFacV1    = new TH2F*[fNcentBins];
 170   fh2PtRawSubFacV1     = new TH2F*[fNcentBins];
 171   fh2PtCorrSubFacV1    = new TH2F*[fNcentBins];
 172   fh2NConstSubFacV1    = new TH2F*[fNcentBins];
 173   fh2PtTrueSubFacV2    = new TH2F*[fNcentBins];
 174   fh2PtRawSubFacV2     = new TH2F*[fNcentBins];
 175   fh2PtCorrSubFacV2    = new TH2F*[fNcentBins];
 176   fh2NConstSubFacV2    = new TH2F*[fNcentBins];
 177
 178   for (Int_t i = 0; i < fNcentBins; i++) {
 179     fh2MSubMatch[i]         = 0;
 180     fh2MSubPtRawAll[i]      = 0;
 181     fh3MSubPtRawDRMatch[i]  = 0;
 182     fh3MSubPtTrueDR[i]      = 0;
 183     fh3MTruePtTrueDR[i]     = 0;
 184     fh3PtTrueDeltaMDR[i]    = 0;
 185     fh3PtTrueDeltaMRelDR[i] = 0;
 186     fhnMassResponse[i]      = 0;
 187     fh2PtTrueSubFacV1[i]    = 0;
 188     fh2PtRawSubFacV1[i]     = 0;
 189     fh2PtCorrSubFacV1[i]    = 0;
 190     fh2NConstSubFacV1[i]    = 0;
 191     fh2PtTrueSubFacV2[i]    = 0;
 192     fh2PtRawSubFacV2[i]     = 0;
 193     fh2PtCorrSubFacV2[i]    = 0;
 194     fh2NConstSubFacV2[i]    = 0;
 195   }
 196
 197   SetMakeGeneralHistograms(kTRUE);
 198   if(fCreateTree) DefineOutput(2, TTree::Class());
 199 }
 200
 201 //________________________________________________________________________
 202 AliAnalysisTaskJetShapeDeriv::~AliAnalysisTaskJetShapeDeriv()
 203 {
 204   // Destructor.
 205 }
 206
 207 //________________________________________________________________________
 208 void AliAnalysisTaskJetShapeDeriv::UserCreateOutputObjects()
 209 {
 210   // Create user output.
 211
 212   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 213
 214   Bool_t oldStatus = TH1::AddDirectoryStatus();
 215   TH1::AddDirectory(kFALSE);
 216
 217   const Int_t nBinsPt  = 200;
 218   const Double_t minPt = -50.;
 219   const Double_t maxPt = 150.;
 220
 221   const Int_t nBinsM  = 100;
 222   const Double_t minM = -20.;
 223   const Double_t maxM = 80.;
 224
 225   // const Int_t nBinsMT  = 50;
 226   // const Double_t minMT = 0.;
 227   // const Double_t maxMT = 50.;
 228
 229   const Int_t nBinsDRToLJ  = 20; //distance to leading jet in Pb-Pb only event
 230   const Double_t minDRToLJ = 0.;
 231   const Double_t maxDRToLJ = 1.;
 232
 233   const Int_t nBinsV1  = 60;
 234   const Double_t minV1 = -60.;
 235   const Double_t maxV1 = 0.;
 236
 237   const Int_t nBinsV2  = 60;
 238   const Double_t minV2 = -30.;
 239   const Double_t maxV2 = 0.;
 240
 241   //Binning for THnSparse
 242   const Int_t nBinsSparse0 = 5;
 243   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
 244   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt, minDRToLJ};
 245   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
 246
 247   TString histName = "";
 248   TString histTitle = "";
 249
 250   for (Int_t i = 0; i < fNcentBins; i++) {
 251     histName = Form("fh2MSubMatch_%d",i);
 252     histTitle = Form("fh2MSubMatch_%d;#it{M}_{sub};match",i);
 253     fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
 254     fOutput->Add(fh2MSubMatch[i]);
 255
 256     histName = Form("fh2MSubPtRawAll_%d",i);
 257     histTitle = Form("fh2MSubPtRawAll_%d;#it{M}_{sub};#it{p}_{T}",i);
 258     fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 259     fOutput->Add(fh2MSubPtRawAll[i]);
 260
 261     histName = Form("fh3MSubPtRawDRMatch_%d",i);
 262     histTitle = Form("fh3MSubPtRawDRMatch_%d;#it{M}_{sub};#it{p}_{T}",i);
 263     fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 264     fOutput->Add(fh3MSubPtRawDRMatch[i]);
 265
 266     histName = Form("fh3MSubPtTrueDR_%d",i);
 267     histTitle = Form("fh3MSubPtTrueDR_%d;#it{M}_{sub};#it{p}_{T}",i);
 268     fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 269     fOutput->Add(fh3MSubPtTrueDR[i]);
 270
 271     histName = Form("fh3MTruePtTrueDR_%d",i);
 272     histTitle = Form("fh3MTruePtTrueDR_%d;#it{M}_{sub};#it{p}_{T}",i);
 273     fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 274     fOutput->Add(fh3MTruePtTrueDR[i]);
 275
 276     histName = Form("fh3PtTrueDeltaMDR_%d",i);
 277     histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#it{M}_{sub}-#it{M}_{true}",i);
 278     fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 279     fOutput->Add(fh3PtTrueDeltaMDR[i]);
 280
 281     histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
 282     histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};(#it{M}_{sub}-#it{M}_{true})/#it{M}_{true}",i);
 283     fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
 284     fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
 285
 286     histName = Form("fhnMassResponse_%d",i);
 287     histTitle = Form("fhnMassResponse_%d;#it{M}_{sub};#it{M}_{true};#it{p}_{T,sub};#it{p}_{T,true};#it{M}_{sub}^{tagged}",i);
 288     fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 289     fOutput->Add(fhnMassResponse[i]);
 290
 291     //derivative histograms
 292     histName = Form("fh2PtTrueSubFacV1_%d",i);
 293     histTitle = Form("fh2PtTrueSubFacV1_%d;#it{p}_{T,true};-(#rho+#rho_{m})V_{1}",i);
 294     fh2PtTrueSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
 295
 296     histName = Form("fh2PtRawSubFacV1_%d",i);
 297     histTitle = Form("fh2PtRawSubFacV1_%d;#it{p}_{T,raw};-(#rho+#rho_{m})V_{1}",i);
 298     fh2PtRawSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
 299
 300     histName = Form("fh2PtCorrSubFacV1_%d",i);
 301     histTitle = Form("fh2PtCorrSubFacV1_%d;#it{p}_{T,corr};-(#rho+#rho_{m})V_{1}",i);
 302     fh2PtCorrSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
 303
 304     histName = Form("fh2NConstSubFacV1_%d",i);
 305     histTitle = Form("fh2NConstSubFacV1_%d;#it{N}_{const};-(#rho+#rho_{m})V_{1}",i);
 306     fh2NConstSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,200.);
 307
 308     histName = Form("fh2PtTrueSubFacV2_%d",i);
 309     histTitle = Form("fh2PtTrueSubFacV2_%d;#it{p}_{T,true};0.5(#rho+#rho_{m})^{2}V_{2}",i);
 310     fh2PtTrueSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
 311
 312     histName = Form("fh2PtRawSubFacV2_%d",i);
 313     histTitle = Form("fh2PtRawSubFacV2_%d;#it{p}_{T,raw};0.5(#rho+#rho_{m})^{2}V_{2}",i);
 314     fh2PtRawSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
 315
 316     histName = Form("fh2PtCorrSubFacV2_%d",i);
 317     histTitle = Form("fh2PtCorrSubFacV2_%d;#it{p}_{T,corr};0.5(#rho+#rho_{m})^{2}V_{2}",i);
 318     fh2PtCorrSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
 319
 320     histName = Form("fh2NConstSubFacV2_%d",i);
 321     histTitle = Form("fh2NConstSubFacV2_%d;#it{N}_{const};0.5(#rho+#rho_{m})^{2}V_{2}",i);
 322     fh2NConstSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,200.);
 323
 324   }
 325
 326   // =========== Switch on Sumw2 for all histos ===========
 327   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 328     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 329     if (h1){
 330       h1->Sumw2();
 331       continue;
 332     }
 333     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 334     if(hn)hn->Sumw2();
 335   }
 336
 337   TH1::AddDirectory(oldStatus);
 338
 339   // Create a tree.
 340   if(fCreateTree) {
 341     fTreeJetBkg = new TTree("fTreeJetBkg", "fTreeJetBkg");
 342     fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
 343     fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
 344     fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
 345     fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
 346     fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
 347     fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
 348     fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
 349     fTreeJetBkg->Branch("fNConst",&fNConst,"fNConst/I");
 350     fTreeJetBkg->Branch("fM1st",&fM1st,"fM1st/F");
 351     fTreeJetBkg->Branch("fM2nd",&fM2nd,"fM2nd/F");
 352     fTreeJetBkg->Branch("fDeriv1st",&fDeriv1st,"fDeriv1st/F");
 353     fTreeJetBkg->Branch("fDeriv2nd",&fDeriv2nd,"fDeriv2nd/F");
 354     fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
 355   }
 356
 357   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 358   if(fCreateTree) PostData(2, fTreeJetBkg);
 359 }
 360
 361 //________________________________________________________________________
 362 Bool_t AliAnalysisTaskJetShapeDeriv::Run()
 363 {
 364   // Run analysis code here, if needed. It will be executed before FillHistograms().
 365
 366   return kTRUE;
 367 }
 368
 369 //________________________________________________________________________
 370 Bool_t AliAnalysisTaskJetShapeDeriv::FillHistograms()
 371 {
 372   // Fill histograms.
 373
 374   AliEmcalJet* jet1  = NULL; //AA jet
 375   AliEmcalJet *jet2  = NULL; //Embedded Pythia jet
 376   //  AliEmcalJet *jet1T = NULL; //tagged AA jet
 377   //  AliEmcalJet *jet2T = NULL; //tagged Pythia jet
 378   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 379   fRho  = (Float_t)jetCont->GetRhoVal();
 380   fRhoM = (Float_t)jetCont->GetRhoMassVal();
 381
 382   //Get leading jet in Pb-Pb event without embedded objects
 383   AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
 384   AliEmcalJet *jetL = NULL;
 385   if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
 386
 387   if(jetCont) {
 388     jetCont->ResetCurrentID();
 389     while((jet1 = jetCont->GetNextAcceptJet())) {
 390       jet2 = NULL;
 391       //      jet1T = NULL;
 392       //   jet2T = NULL;
 393       //      if(jet1->GetTagStatus()>0) jet1T = jet1->GetTaggedJet();
 394       //Fill histograms for all AA jets
 395       fh2MSubPtRawAll[fCentBin]->Fill(jet1->GetSecondOrderSubtracted(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
 396       fh2PtRawSubFacV1[fCentBin]->Fill(jet1->Pt(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
 397       fh2PtCorrSubFacV1[fCentBin]->Fill(jet1->Pt()-fRho*jet1->Area(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
 398       fh2NConstSubFacV1[fCentBin]->Fill(jet1->GetNumberOfTracks(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
 399       fh2PtRawSubFacV2[fCentBin]->Fill(jet1->Pt(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
 400       fh2PtCorrSubFacV2[fCentBin]->Fill(jet1->Pt()-fRho*jet1->Area(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
 401       fh2NConstSubFacV2[fCentBin]->Fill(jet1->GetNumberOfTracks(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
 402
 403       Double_t fraction = 0.;
 404       fMatch = 0;
 405       fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
 406       if(fSingleTrackEmb) {
 407         AliVParticle *vp = GetEmbeddedConstituent(jet1);
 408         if(vp) {
 409           fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 410           fMatch = 1;
 411         }
 412       } else {
 413         jet2 = jet1->ClosestJet();
 414         fraction = jetCont->GetFractionSharedPt(jet1);
 415         fMatch = 1;
 416         if(fMinFractionShared>0.) {
 417           if(fraction>fMinFractionShared) {
 418             fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
 419             fMatch = 1;
 420           } else
 421             fMatch = 0;
 422         }
 423       }
 424
 425       //      if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
 426
 427       //Fill histograms for matched jets
 428       fh2MSubMatch[fCentBin]->Fill(jet1->GetSecondOrderSubtracted(),fMatch);
 429       if(fMatch==1) {
 430         Double_t drToLJ = -1.;
 431         if(jetL) drToLJ = jet1->DeltaR(jetL);
 432         fh3MSubPtRawDRMatch[fCentBin]->Fill(jet1->GetSecondOrderSubtracted(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area(),drToLJ);
 433         if(jet2) {
 434           fh3MSubPtTrueDR[fCentBin]->Fill(jet1->GetSecondOrderSubtracted(),jet2->Pt(),drToLJ);
 435           fh3MTruePtTrueDR[fCentBin]->Fill(jet2->M(),jet2->Pt(),drToLJ);
 436           fh3PtTrueDeltaMDR[fCentBin]->Fill(jet2->Pt(),jet1->GetSecondOrderSubtracted()-jet2->M(),drToLJ);
 437           if(jet2->M()>0.) fh3PtTrueDeltaMRelDR[fCentBin]->Fill(jet2->Pt(),(jet1->GetSecondOrderSubtracted()-jet2->M())/jet2->M(),drToLJ);
 438           // Double_t mJet1Tagged = -1.;
 439           // if(jet1T) mJet1Tagged = jet1T->M();
 440           Double_t var[5] = {jet1->GetSecondOrderSubtracted(),jet2->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area(),jet2->Pt(),drToLJ};
 441           fhnMassResponse[fCentBin]->Fill(var);
 442         }
 443       }
 444
 445       if(fCreateTree) {
 446         fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
 447         fArea = (Float_t)jet1->Area();
 448         fAreaPhi = (Float_t)jet1->AreaPhi();
 449         fAreaEta = (Float_t)jet1->AreaEta();
 450         fNConst = (Int_t)jet1->GetNumberOfTracks();
 451         fM1st   = (Float_t)jet1->GetFirstOrderSubtracted();
 452         fM2nd   = (Float_t)jet1->GetSecondOrderSubtracted();
 453         fDeriv1st = (Float_t)jet1->GetFirstDerivative();
 454         fDeriv2nd = (Float_t)jet1->GetSecondDerivative();
 455         fTreeJetBkg->Fill();
 456       }
 457     }
 458   }
 459   return kTRUE;
 460 }
 461
 462 //________________________________________________________________________
 463 AliVParticle* AliAnalysisTaskJetShapeDeriv::GetEmbeddedConstituent(AliEmcalJet *jet) {
 464
 465   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 466   AliVParticle *vp = 0x0;
 467   AliVParticle *vpe = 0x0; //embedded particle
 468   Int_t nc = 0;
 469   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 470     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
 471     if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
 472     if(!vpe) vpe = vp;
 473     else if(vp->Pt()>vpe->Pt()) vpe = vp;
 474     nc++;
 475   }
 476   AliDebug(11,Form("Found %d embedded particles",nc));
 477   return vpe;
 478 }
 479
 480
 481 //________________________________________________________________________
 482 Bool_t AliAnalysisTaskJetShapeDeriv::RetrieveEventObjects() {
 483   //
 484   // retrieve event objects
 485   //
 486
 487   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 488     return kFALSE;
 489
 490   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 491   jetCont->LoadRhoMass(InputEvent());
 492
 493   return kTRUE;
 494 }
 495
 496 //_______________________________________________________________________
 497 void AliAnalysisTaskJetShapeDeriv::Terminate(Option_t *)
 498 {
 499   // Called once at the end of the analysis.
 500 }
 501