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