PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMassResponseDet.cxx

   1 //
   2 // Detector response jet mass analysis task.
   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 <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TProfile.h>
  14 #include <TChain.h>
  15 #include <TSystem.h>
  16 #include <TFile.h>
  17 #include <TKey.h>
  18
  19 #include "AliVCluster.h"
  20 #include "AliVTrack.h"
  21 #include "AliEmcalJet.h"
  22 #include "AliRhoParameter.h"
  23 #include "AliLog.h"
  24 #include "AliEmcalParticle.h"
  25 #include "AliMCEvent.h"
  26 #include "AliGenPythiaEventHeader.h"
  27 #include "AliAODMCHeader.h"
  28 #include "AliMCEvent.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliJetContainer.h"
  31
  32 #include "AliAODEvent.h"
  33
  34 #include "AliAnalysisTaskJetMassResponseDet.h"
  35
  36 ClassImp(AliAnalysisTaskJetMassResponseDet)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskJetMassResponseDet::AliAnalysisTaskJetMassResponseDet() :
  40   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetMassResponseDet", kTRUE),
  41   fContainerPart(0),
  42   fContainerDet(0),
  43   fJetMassType(kRaw),
  44   fh2PtVsMassJetPartAll(0),
  45   fh2PtVsMassJetPartMatch(0),
  46   fh2PtVsMassJetPartTagged(0),
  47   fh2PtVsMassJetPartTaggedMatch(0),
  48   fh2PtVsMassJetDetAll(0),
  49   fh2PtVsMassJetDetTagged(0),
  50   fh2EtaPhiMatchedDet(0),
  51   fh2EtaPhiMatchedPart(0),
  52   fhnMassResponse(0),
  53   fh1AreaPartAll(0),
  54   fh1AreaDetAll(0)
  55 {
  56   // Default constructor.
  57
  58   SetMakeGeneralHistograms(kTRUE);
  59 }
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskJetMassResponseDet::AliAnalysisTaskJetMassResponseDet(const char *name) :
  63   AliAnalysisTaskEmcalJet(name, kTRUE),
  64   fContainerPart(0),
  65   fContainerDet(0),
  66   fJetMassType(kRaw),
  67   fh2PtVsMassJetPartAll(0),
  68   fh2PtVsMassJetPartMatch(0),
  69   fh2PtVsMassJetPartTagged(0),
  70   fh2PtVsMassJetPartTaggedMatch(0),
  71   fh2PtVsMassJetDetAll(0),
  72   fh2PtVsMassJetDetTagged(0),
  73   fh2EtaPhiMatchedDet(0),
  74   fh2EtaPhiMatchedPart(0),
  75   fhnMassResponse(0),
  76   fh1AreaPartAll(0),
  77   fh1AreaDetAll(0)
  78 {
  79   // Standard constructor.
  80
  81   SetMakeGeneralHistograms(kTRUE);
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskJetMassResponseDet::~AliAnalysisTaskJetMassResponseDet()
  86 {
  87   // Destructor.
  88 }
  89
  90 //________________________________________________________________________
  91 void AliAnalysisTaskJetMassResponseDet::UserCreateOutputObjects()
  92 {
  93   // Create user output.
  94
  95   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  96
  97   Bool_t oldStatus = TH1::AddDirectoryStatus();
  98   TH1::AddDirectory(kFALSE);
  99
 100   const Int_t nBinsPt  = 200;
 101   const Double_t minPt = 0.;
 102   const Double_t maxPt = 200.;
 103
 104   const Int_t nBinsM  = 200;
 105   const Double_t minM = 0.;
 106   const Double_t maxM = 50.;
 107
 108   const Int_t nBinsConstEff  = 40;
 109   const Double_t minConstEff = 0.;
 110   const Double_t maxConstEff = 2.;
 111
 112   // const Int_t nBinsConst = 26;
 113   // const Double_t minConst = -5.5;
 114   // const Double_t maxConst = 20.5;
 115
 116   //Binning for THnSparse
 117   const Int_t nBinsSparse0 = 5;
 118   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsConstEff};
 119   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt, minConstEff};
 120   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt, maxConstEff};
 121
 122   //Create histograms
 123   TString histName = "";
 124   TString histTitle = "";
 125
 126   histName = "fh2PtVsMassJetPartAll";
 127   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 128   fh2PtVsMassJetPartAll = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 129   fOutput->Add(fh2PtVsMassJetPartAll);
 130
 131   histName = "fh2PtVsMassJetPartMatch";
 132   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 133   fh2PtVsMassJetPartMatch = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 134   fOutput->Add(fh2PtVsMassJetPartMatch);
 135
 136   histName = "fh2PtVsMassJetPartTagged";
 137   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 138   fh2PtVsMassJetPartTagged = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 139   fOutput->Add(fh2PtVsMassJetPartTagged);
 140
 141   histName = "fh2PtVsMassJetPartTaggedMatch";
 142   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 143   fh2PtVsMassJetPartTaggedMatch = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 144   fOutput->Add(fh2PtVsMassJetPartTaggedMatch);
 145
 146   histName = "fh2PtVsMassJetDetAll";
 147   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 148   fh2PtVsMassJetDetAll = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 149   fOutput->Add(fh2PtVsMassJetDetAll);
 150
 151   histName = "fh2PtVsMassJetDetTagged";
 152   histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
 153   fh2PtVsMassJetDetTagged = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 154   fOutput->Add(fh2PtVsMassJetDetTagged);
 155
 156   histName = "fh2EtaPhiMatchedDet";
 157   histTitle = TString::Format("%s;#eta;#varphi",histName.Data());
 158   fh2EtaPhiMatchedDet = new TH2F(histName.Data(),histTitle.Data(),100,-1.,1.,72,0.,TMath::TwoPi());
 159   fOutput->Add(fh2EtaPhiMatchedDet);
 160
 161   histName = "fh2EtaPhiMatchedPart";
 162   histTitle = TString::Format("%s;#eta;#varphi",histName.Data());
 163   fh2EtaPhiMatchedPart = new TH2F(histName.Data(),histTitle.Data(),100,-1.,1.,72,0.,TMath::TwoPi());
 164   fOutput->Add(fh2EtaPhiMatchedPart);
 165
 166   histName = "fhnMassResponse";
 167   histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part};#it{N}_{const}^{det}/#it{N}_{const}^{part}",histName.Data());
 168   fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 169   fOutput->Add(fhnMassResponse);
 170
 171   fh1AreaPartAll = new TH1D("fh1AreaPartAll","fh1AreaPartAll",100.,0.,1.);
 172   fOutput->Add(fh1AreaPartAll);
 173
 174   fh1AreaDetAll = new TH1D("fh1AreaDetAll","fh1AreaDetAll",100.,0.,1.);
 175   fOutput->Add(fh1AreaDetAll);
 176
 177
 178   // =========== Switch on Sumw2 for all histos ===========
 179   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 180     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 181     if (h1){
 182       h1->Sumw2();
 183       continue;
 184     }
 185     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 186     if(hn)hn->Sumw2();
 187   }
 188
 189   TH1::AddDirectory(oldStatus);
 190
 191   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 192 }
 193
 194 //________________________________________________________________________
 195 Bool_t AliAnalysisTaskJetMassResponseDet::Run()
 196 {
 197   // Run analysis code here, if needed. It will be executed before FillHistograms().
 198   return kTRUE;
 199 }
 200
 201 //________________________________________________________________________
 202 Bool_t AliAnalysisTaskJetMassResponseDet::FillHistograms()
 203 {
 204   // Fill histograms.
 205
 206   AliJetContainer *cPart = GetJetContainer(fContainerPart);
 207   AliJetContainer *cDet = GetJetContainer(fContainerDet);
 208   AliEmcalJet* jPart = NULL;
 209   AliEmcalJet* jDet = NULL;
 210
 211   //loop on particle level jets
 212   Int_t nAccPart = 0;
 213   if(cPart) {
 214     cPart->ResetCurrentID();
 215     while((jPart = cPart->GetNextAcceptJet())) {
 216       fh2PtVsMassJetPartAll->Fill(jPart->Pt(),jPart->M());
 217       fh1AreaPartAll->Fill(jPart->Area());
 218       nAccPart++;
 219       jDet = jPart->ClosestJet();
 220       if(jDet) fh2PtVsMassJetPartMatch->Fill(jPart->Pt(),jPart->M());
 221       if(jPart->GetTagStatus()<1 || !jPart->GetTaggedJet())
 222         continue;
 223       fh2PtVsMassJetPartTagged->Fill(jPart->Pt(),jPart->M());
 224       if(jDet) fh2PtVsMassJetPartTaggedMatch->Fill(jPart->Pt(),jPart->M());
 225     }
 226   }
 227
 228   //loop on detector level jets
 229   Int_t nAccDet = 0;
 230   if(cDet) {
 231     cDet->ResetCurrentID();
 232     while((jDet = cDet->GetNextAcceptJet())) {
 233       Double_t mjet = GetJetMass(jDet);
 234       fh2PtVsMassJetDetAll->Fill(jDet->Pt(),mjet);
 235       fh1AreaDetAll->Fill(jDet->Area());
 236       nAccDet++;
 237       if(jDet->GetTagStatus()>=1 && jDet->GetTaggedJet())
 238          fh2PtVsMassJetDetTagged->Fill(jDet->Pt(),mjet);
 239
 240        //fill detector response
 241        jPart = jDet->ClosestJet();
 242        if(jPart) {
 243          fh2EtaPhiMatchedDet->Fill(jDet->Eta(),jDet->Phi());
 244          fh2EtaPhiMatchedPart->Fill(jPart->Eta(),jPart->Phi());
 245
 246          Int_t nConstPart = jPart->GetNumberOfConstituents();
 247          Int_t nConstDet = jDet->GetNumberOfConstituents();
 248          Double_t eff = -1.;
 249          if(nConstPart>0) eff = (Double_t)nConstDet/((Double_t)nConstPart);
 250          Double_t var[5] = {GetJetMass(jDet),jPart->M(),jDet->Pt(),jPart->Pt(),eff};
 251          fhnMassResponse->Fill(var);
 252        }
 253     }
 254   }
 255
 256   return kTRUE;
 257 }
 258
 259 //________________________________________________________________________
 260 Double_t AliAnalysisTaskJetMassResponseDet::GetJetMass(AliEmcalJet *jet) {
 261   //calc subtracted jet mass
 262   if(fJetMassType==kRaw)
 263     return jet->M();
 264   else if(fJetMassType==kDeriv)
 265     return jet->GetSecondOrderSubtracted();
 266
 267   return 0;
 268 }
 269
 270 //________________________________________________________________________
 271 Bool_t AliAnalysisTaskJetMassResponseDet::RetrieveEventObjects() {
 272   //
 273   // retrieve event objects
 274   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 275     return kFALSE;
 276
 277   return kTRUE;
 278 }
 279
 280 //_______________________________________________________________________
 281 void AliAnalysisTaskJetMassResponseDet::Terminate(Option_t *)
 282 {
 283   // Called once at the end of the analysis.
 284 }
 285