PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeConst.cxx

   1 //
   2 // Do subtraction for jet shapes using constituents arXiv:1403.3108
   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 "AliAnalysisTaskJetShapeConst.h"
  36
  37 ClassImp(AliAnalysisTaskJetShapeConst)
  38
  39 //________________________________________________________________________
  40 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() :
  41   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
  42   fContainerBase(0),
  43   fContainerSub(1),
  44   fContainerNoEmb(2),
  45   fMinFractionShared(0),
  46   fSingleTrackEmb(kFALSE),
  47   fCreateTree(kFALSE),
  48   fTreeJetBkg(),
  49   fJet1Vec(new TLorentzVector()),
  50   fJet2Vec(new TLorentzVector()),
  51   fJetSubVec(new TLorentzVector()),
  52   fArea(0),
  53   fAreaPhi(0),
  54   fAreaEta(0),
  55   fRho(0),
  56   fRhoM(0),
  57   fNConst(0),
  58   fMatch(0),
  59   fh2MSubMatch(0x0),
  60   fh2MSubPtRawAll(0x0),
  61   fh2MSubPtRawMatch(0x0),
  62   fh2MSubPtTrue(0x0),
  63   fh2MTruePtTrue(0x0),
  64   fh2PtTrueDeltaM(0x0),
  65   fh2PtTrueDeltaMRel(0x0),
  66   fhnMassResponse(0x0)
  67 {
  68   // Default constructor.
  69
  70   fh2MSubMatch        = new TH2F*[fNcentBins];
  71   fh2MSubPtRawAll     = new TH2F*[fNcentBins];
  72   fh2MSubPtRawMatch   = new TH2F*[fNcentBins];
  73   fh2MSubPtTrue       = new TH2F*[fNcentBins];
  74   fh2MTruePtTrue      = new TH2F*[fNcentBins];
  75   fh2PtTrueDeltaM     = new TH2F*[fNcentBins];
  76   fh2PtTrueDeltaMRel  = new TH2F*[fNcentBins];
  77   fhnMassResponse     = new THnSparse*[fNcentBins];
  78
  79   for (Int_t i = 0; i < fNcentBins; i++) {
  80     fh2MSubMatch[i]        = 0;
  81     fh2MSubPtRawAll[i]     = 0;
  82     fh2MSubPtRawMatch[i]   = 0;
  83     fh2MSubPtTrue[i]       = 0;
  84     fh2MTruePtTrue[i]      = 0;
  85     fh2PtTrueDeltaM[i]     = 0;
  86     fh2PtTrueDeltaMRel[i]  = 0;
  87     fhnMassResponse[i]     = 0;
  88   }
  89
  90   SetMakeGeneralHistograms(kTRUE);
  91   if(fCreateTree) DefineOutput(2, TTree::Class());
  92 }
  93
  94 //________________________________________________________________________
  95 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
  96   AliAnalysisTaskEmcalJet(name, kTRUE),
  97   fContainerBase(0),
  98   fContainerSub(1),
  99   fContainerNoEmb(2),
 100   fMinFractionShared(0),
 101   fSingleTrackEmb(kFALSE),
 102   fCreateTree(kFALSE),
 103   fTreeJetBkg(0),
 104   fJet1Vec(new TLorentzVector()),
 105   fJet2Vec(new TLorentzVector()),
 106   fJetSubVec(new TLorentzVector()),
 107   fArea(0),
 108   fAreaPhi(0),
 109   fAreaEta(0),
 110   fRho(0),
 111   fRhoM(0),
 112   fNConst(0),
 113   fMatch(0),
 114   fh2MSubMatch(0x0),
 115   fh2MSubPtRawAll(0x0),
 116   fh2MSubPtRawMatch(0x0),
 117   fh2MSubPtTrue(0x0),
 118   fh2MTruePtTrue(0x0),
 119   fh2PtTrueDeltaM(0x0),
 120   fh2PtTrueDeltaMRel(0x0),
 121   fhnMassResponse(0x0)
 122 {
 123   // Standard constructor.
 124
 125   fh2MSubMatch        = new TH2F*[fNcentBins];
 126   fh2MSubPtRawAll     = new TH2F*[fNcentBins];
 127   fh2MSubPtRawMatch   = new TH2F*[fNcentBins];
 128   fh2MSubPtTrue       = new TH2F*[fNcentBins];
 129   fh2MTruePtTrue      = new TH2F*[fNcentBins];
 130   fh2PtTrueDeltaM     = new TH2F*[fNcentBins];
 131   fh2PtTrueDeltaMRel  = new TH2F*[fNcentBins];
 132   fhnMassResponse     = new THnSparse*[fNcentBins];
 133
 134   for (Int_t i = 0; i < fNcentBins; i++) {
 135     fh2MSubMatch[i]        = 0;
 136     fh2MSubPtRawAll[i]     = 0;
 137     fh2MSubPtRawMatch[i]   = 0;
 138     fh2MSubPtTrue[i]       = 0;
 139     fh2MTruePtTrue[i]      = 0;
 140     fh2PtTrueDeltaM[i]     = 0;
 141     fh2PtTrueDeltaMRel[i]  = 0;
 142     fhnMassResponse[i]     = 0;
 143   }
 144
 145   SetMakeGeneralHistograms(kTRUE);
 146   if(fCreateTree) DefineOutput(2, TTree::Class());
 147 }
 148
 149 //________________________________________________________________________
 150 AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
 151 {
 152   // Destructor.
 153 }
 154
 155 //________________________________________________________________________
 156 void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
 157 {
 158   // Create user output.
 159
 160   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 161
 162   Bool_t oldStatus = TH1::AddDirectoryStatus();
 163   TH1::AddDirectory(kFALSE);
 164
 165   const Int_t nBinsPt  = 200;
 166   const Double_t minPt = -50.;
 167   const Double_t maxPt = 150.;
 168
 169   const Int_t nBinsM  = 100;
 170   const Double_t minM = -20.;
 171   const Double_t maxM = 80.;
 172
 173   // const Int_t nBinsMT  = 50;
 174   // const Double_t minMT = 0.;
 175   // const Double_t maxMT = 50.;
 176
 177   const Int_t nBinsDRToLJ  = 20; //distance to leading jet in Pb-Pb only event
 178   const Double_t minDRToLJ = 0.;
 179   const Double_t maxDRToLJ = 1.;
 180
 181   //Binning for THnSparse
 182   const Int_t nBinsSparse0 = 5;
 183   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
 184   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt, minDRToLJ};
 185   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
 186
 187   TString histName = "";
 188   TString histTitle = "";
 189   for (Int_t i = 0; i < fNcentBins; i++) {
 190     histName = Form("fh2MSubMatch_%d",i);
 191     histTitle = Form("fh2MSubMatch_%d;#it{M}_{sub};match",i);
 192     fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
 193     fOutput->Add(fh2MSubMatch[i]);
 194
 195     histName = Form("fh2MSubPtRawAll_%d",i);
 196     histTitle = Form("fh2MSubPtRawAll_%d;#it{M}_{sub};#it{p}_{T}",i);
 197     fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 198     fOutput->Add(fh2MSubPtRawAll[i]);
 199
 200     histName = Form("fh2MSubPtRawMatch_%d",i);
 201     histTitle = Form("fh2MSubPtRawMatch_%d;#it{M}_{sub};#it{p}_{T}",i);
 202     fh2MSubPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 203     fOutput->Add(fh2MSubPtRawMatch[i]);
 204
 205     histName = Form("fh2MSubPtTrue_%d",i);
 206     histTitle = Form("fh2MSubPtTrue_%d;#it{M}_{sub};#it{p}_{T}",i);
 207     fh2MSubPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 208     fOutput->Add(fh2MSubPtTrue[i]);
 209
 210     histName = Form("fh2MTruePtTrue_%d",i);
 211     histTitle = Form("fh2MTruePtTrue_%d;#it{M}_{sub};#it{p}_{T}",i);
 212     fh2MTruePtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
 213     fOutput->Add(fh2MTruePtTrue[i]);
 214
 215     histName = Form("fh2PtTrueDeltaM_%d",i);
 216     histTitle = Form("fh2PtTrueDeltaM_%d;#it{p}_{T,true};#it{M}_{sub}-#it{M}_{true}",i);
 217     fh2PtTrueDeltaM[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.);
 218     fOutput->Add(fh2PtTrueDeltaM[i]);
 219
 220     histName = Form("fh2PtTrueDeltaMRel_%d",i);
 221     histTitle = Form("fh2PtTrueDeltaMRel_%d;#it{p}_{T,true};(#it{M}_{sub}-#it{M}_{true})/#it{M}_{true}",i);
 222     fh2PtTrueDeltaMRel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.);
 223     fOutput->Add(fh2PtTrueDeltaMRel[i]);
 224
 225     histName = Form("fhnMassResponse_%d",i);
 226     histTitle = Form("fhnMassResponse_%d;#it{M}_{sub};#it{M}_{true};#it{p}_{T,sub};#it{p}_{T,true};#it{M}_{sub}^{tagged}",i);
 227     fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 228     fOutput->Add(fhnMassResponse[i]);
 229   }
 230
 231   // =========== Switch on Sumw2 for all histos ===========
 232   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 233     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 234     if (h1){
 235       h1->Sumw2();
 236       continue;
 237     }
 238     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 239     if(hn)hn->Sumw2();
 240   }
 241
 242   TH1::AddDirectory(oldStatus);
 243
 244   // Create a tree.
 245   if(fCreateTree) {
 246     fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
 247     fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
 248     fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
 249     fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
 250     fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
 251     fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
 252     fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
 253     fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
 254     fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
 255     fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
 256   }
 257
 258   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 259   if(fCreateTree) PostData(2, fTreeJetBkg);
 260 }
 261
 262 //________________________________________________________________________
 263 Bool_t AliAnalysisTaskJetShapeConst::Run()
 264 {
 265   // Run analysis code here, if needed. It will be executed before FillHistograms().
 266
 267   return kTRUE;
 268 }
 269
 270 //________________________________________________________________________
 271 Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
 272 {
 273   // Fill histograms.
 274
 275   AliEmcalJet* jet1  = NULL; //AA jet
 276   AliEmcalJet *jet2  = NULL; //Embedded Pythia jet
 277   //  AliEmcalJet *jet1T = NULL; //tagged AA jet
 278   //  AliEmcalJet *jet2T = NULL; //tagged Pythia jet
 279   AliEmcalJet *jetS = NULL;  //subtracted jet
 280   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 281   AliJetContainer *jetContS = GetJetContainer(fContainerSub);
 282
 283   //Get leading jet in Pb-Pb event without embedded objects
 284   AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
 285   AliEmcalJet *jetL = NULL;
 286   if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
 287
 288   AliDebug(11,Form("NJets  Incl: %d  Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
 289   if(jetCont) {
 290     jetCont->ResetCurrentID();
 291     while((jet1 = jetCont->GetNextAcceptJet())) {
 292       jet2  = NULL;
 293       //      jet1T = NULL;
 294       //   jet2T = NULL;
 295       jetS  = NULL;
 296       //      if(jet1->GetTagStatus()>0) jet1T = jet1->GetTaggedJet();
 297
 298       //Get constituent subtacted version of jet
 299       Int_t ifound = 0;
 300       Int_t ilab = -1;
 301       for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
 302         //if(ifound==1) continue;
 303         jetS = jetContS->GetJet(i);
 304         if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
 305           ifound++;
 306           if(ifound==1) ilab = i;
 307         }
 308       }
 309       if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
 310       if(ifound==0) jetS = 0x0;
 311       else jetS = jetContS->GetJet(ilab);
 312
 313       //Fill histograms for all AA jets
 314       fh2MSubPtRawAll[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
 315
 316       Double_t fraction = 0.;
 317       fMatch = 0;
 318       fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
 319       if(fSingleTrackEmb) {
 320         AliVParticle *vp = GetEmbeddedConstituent(jet1);
 321         if(vp) {
 322           fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
 323           fMatch = 1;
 324         }
 325       } else {
 326         jet2 = jet1->ClosestJet();
 327         fraction = jetCont->GetFractionSharedPt(jet1);
 328         fMatch = 1;
 329         if(fMinFractionShared>0.) {
 330           if(fraction>fMinFractionShared) {
 331             fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
 332             fMatch = 1;
 333           } else
 334             fMatch = 0;
 335         }
 336       }
 337
 338       //      if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
 339
 340       //Fill histograms for matched jets
 341       fh2MSubMatch[fCentBin]->Fill(jetS->M(),fMatch);
 342       if(fMatch==1) {
 343         fh2MSubPtRawMatch[fCentBin]->Fill(jetS->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area());
 344         if(jet2) {
 345           fh2MSubPtTrue[fCentBin]->Fill(jetS->M(),jet2->Pt());
 346           fh2MTruePtTrue[fCentBin]->Fill(jet2->M(),jet2->Pt());
 347           fh2PtTrueDeltaM[fCentBin]->Fill(jet2->Pt(),jetS->M()-jet2->M());
 348           if(jet2->M()>0.) fh2PtTrueDeltaMRel[fCentBin]->Fill(jet2->Pt(),(jetS->M()-jet2->M())/jet2->M());
 349           //      Double_t mJet1Tagged = -1.;
 350           //      if(jet1T) mJet1Tagged = jet1T->M();
 351           Double_t drToLJ = 0.;
 352           if(jetL) drToLJ = jet1->DeltaR(jetL);
 353           Double_t var[5] = {jetS->M(),jet2->M(),jet1->Pt()-jetCont->GetRhoVal()*jet1->Area(),jet2->Pt(),drToLJ};
 354           fhnMassResponse[fCentBin]->Fill(var);
 355         }
 356       }
 357
 358       if(fCreateTree) {
 359         fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
 360         if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
 361         else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
 362         fArea = (Float_t)jet1->Area();
 363         fAreaPhi = (Float_t)jet1->AreaPhi();
 364         fAreaEta = (Float_t)jet1->AreaEta();
 365         fRho  = (Float_t)jetCont->GetRhoVal();
 366         fRhoM = (Float_t)jetCont->GetRhoMassVal();
 367         fNConst = (Int_t)jet1->GetNumberOfTracks();
 368         fTreeJetBkg->Fill();
 369       }
 370     } //jet1 loop
 371   }//jetCont
 372
 373
 374   return kTRUE;
 375 }
 376
 377 //________________________________________________________________________
 378 AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
 379
 380   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 381   AliVParticle *vp = 0x0;
 382   AliVParticle *vpe = 0x0; //embedded particle
 383   Int_t nc = 0;
 384   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
 385     vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
 386     if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
 387     if(!vpe) vpe = vp;
 388     else if(vp->Pt()>vpe->Pt()) vpe = vp;
 389     nc++;
 390   }
 391   AliDebug(11,Form("Found %d embedded particles",nc));
 392   return vpe;
 393 }
 394
 395
 396 //________________________________________________________________________
 397 Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
 398   //
 399   // retrieve event objects
 400   //
 401
 402   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 403     return kFALSE;
 404
 405   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 406   jetCont->LoadRhoMass(InputEvent());
 407
 408   return kTRUE;
 409 }
 410
 411 //_______________________________________________________________________
 412 void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
 413 {
 414   // Called once at the end of the analysis.
 415 }
 416