PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetAna.cxx

   1 //
   2 // Dijet 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 "AliAnalysisUtils.h"
  25 #include "AliEmcalParticle.h"
  26 #include "AliMCEvent.h"
  27 #include "AliGenPythiaEventHeader.h"
  28 #include "AliAODMCHeader.h"
  29 #include "AliMCEvent.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliJetContainer.h"
  32
  33 #include "AliAnalysisTaskEmcalDiJetAna.h"
  34
  35 ClassImp(AliAnalysisTaskEmcalDiJetAna)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna() :
  39   AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetAna"),
  40   fDoMatchFullCharged(kTRUE),
  41   fh2CentRhoCh(0),
  42   fh2CentRhoScaled(0),
  43   fh3PtEtaPhiJetFull(0),
  44   fh3PtEtaPhiJetCharged(0),
  45   fhnDiJetVarsFull(0),
  46   fhnDiJetVarsCh(0),
  47   fhnDiJetVarsFullCharged(0),
  48   fhnMatchingFullCharged(0),
  49   fh3JetPtFullFractionDR(0)
  50 {
  51   // Default constructor.
  52
  53   SetMakeGeneralHistograms(kTRUE);
  54 }
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna(const char *name) :
  58   AliAnalysisTaskEmcalDiJetBase(name),
  59   fDoMatchFullCharged(kTRUE),
  60   fh2CentRhoCh(0),
  61   fh2CentRhoScaled(0),
  62   fh3PtEtaPhiJetFull(0),
  63   fh3PtEtaPhiJetCharged(0),
  64   fhnDiJetVarsFull(0),
  65   fhnDiJetVarsCh(0),
  66   fhnDiJetVarsFullCharged(0),
  67   fhnMatchingFullCharged(0),
  68   fh3JetPtFullFractionDR(0)
  69 {
  70   // Standard constructor.
  71
  72   SetMakeGeneralHistograms(kTRUE);
  73 }
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskEmcalDiJetAna::~AliAnalysisTaskEmcalDiJetAna()
  77 {
  78   // Destructor.
  79 }
  80
  81 //________________________________________________________________________
  82 Bool_t AliAnalysisTaskEmcalDiJetAna::RetrieveEventObjects() {
  83   //
  84   // retrieve event objects
  85   //
  86
  87   if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
  88     return kFALSE;
  89
  90   return kTRUE;
  91 }
  92
  93 //________________________________________________________________________
  94 void AliAnalysisTaskEmcalDiJetAna::UserCreateOutputObjects()
  95 {
  96   // Create user output.
  97
  98   InitOnce();
  99
 100   AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
 101
 102   Bool_t oldStatus = TH1::AddDirectoryStatus();
 103   TH1::AddDirectory(kFALSE);
 104
 105   const Int_t nBinsCent = 100.;
 106   Double_t minCent = 0.;
 107   Double_t maxCent = 100.;
 108
 109   const Int_t nBinsRho = 200;
 110   Double_t minRho = 0.;
 111   Double_t maxRho = 20.;
 112   fh2CentRhoCh = new TH2F("fh2CentRhoCh","fh2CentRhoCh;centrality;#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
 113   fOutput->Add(fh2CentRhoCh);
 114   fh2CentRhoScaled = new TH2F("fh2CentRhoScaled","fh2CentRhoScaled;centrality;s_{EMC}#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
 115   fOutput->Add(fh2CentRhoScaled);
 116
 117   const Int_t nBinsPt = 120;
 118   Double_t minPt = -20.;
 119   Double_t maxPt = 100.;
 120   const Int_t nBinsEta = 40;
 121   Double_t minEta = -1.;
 122   Double_t maxEta = 1.;
 123   const Int_t nBinsPhi = 18*6;
 124   Double_t minPhi = 0.;
 125   Double_t maxPhi = TMath::TwoPi();
 126
 127   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 128   fOutput->Add(fh3PtEtaPhiJetFull);
 129
 130   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 131   fOutput->Add(fh3PtEtaPhiJetCharged);
 132
 133   const Int_t nBinsSparse0 = 6;
 134   const Int_t nBinsDPhi     = 72;
 135   const Int_t nBinsKt       = 100;
 136   const Int_t nBinsDiJetEta = 40;
 137   const Int_t nBinsCentr    = 10;
 138   const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsDPhi,nBinsKt,nBinsDiJetEta,nBinsCentr};
 139   //pT1, pT2, deltaPhi, kT
 140   const Double_t xmin0[nBinsSparse0]  = {  minPt, minPt, -0.5*TMath::Pi(),-100.,-1.,0.};
 141   const Double_t xmax0[nBinsSparse0]  = {  maxPt, maxPt,  1.5*TMath::Pi(), 100., 1.,100.};
 142
 143   if(fDoChargedCharged) {
 144     fhnDiJetVarsCh = new THnSparseF("fhnDiJetVarsCh",
 145                                     "fhnDiJetVarsCh;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality",
 146                                     nBinsSparse0,nBins0,xmin0,xmax0);
 147     fOutput->Add(fhnDiJetVarsCh);
 148   }
 149
 150   if(fDoFullCharged) {
 151     fhnDiJetVarsFullCharged = new THnSparseF("fhnDiJetVarsFullCharged",
 152                                 "fhnDiJetVarsFullCharged;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality",
 153                                 nBinsSparse0,nBins0,xmin0,xmax0);
 154     fOutput->Add(fhnDiJetVarsFullCharged);
 155   }
 156
 157   const Int_t nBinsSparseMatch = 7;
 158   const Int_t nBinsDPhiMatch = 80;
 159   const Int_t nBinsDEtaMatch = 80;
 160   const Int_t nBinsDR = 20;
 161   const Int_t nBinsFraction = 21;
 162   const Int_t nBinsType = 3;
 163   const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsFraction,nBinsType};
 164   //pTfull, pTch, deltaPhi, deltaEta, deltaR, fraction, jet type (leading,subleading,other)
 165   const Double_t xminMatch[nBinsSparseMatch]  = { minPt, minPt, -0.5,-0.5, 0.,0.  ,0};
 166   const Double_t xmaxMatch[nBinsSparseMatch]  = { maxPt, maxPt,  0.5, 0.5, 0.5,1.05,3};
 167   fhnMatchingFullCharged = new THnSparseF("fhnMatchingFullCharged","fhnMatchingFullCharged;#it{p}_{T,full} (GeV/#it{c});#it{p}_{T,ch} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;f_{ch};type",
 168                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 169   fOutput->Add(fhnMatchingFullCharged);
 170
 171   fh3JetPtFullFractionDR = new TH3F("fh3JetPtFullFractionDR","fh3JetPtFullFractionDR;#it{p}_{T,full} (GeV/#it{c}); #it{f}_{ch};#Delta R",nBinsPt,minPt,maxPt,nBinsFraction,0.,1.05,nBinsDR,0.,1.);
 172   fOutput->Add(fh3JetPtFullFractionDR);
 173
 174   // =========== Switch on Sumw2 for all histos ===========
 175   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 176     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 177     if (h1){
 178       h1->Sumw2();
 179       continue;
 180     }
 181     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 182     if(hn) {
 183       hn->Sumw2();
 184       continue;
 185     }
 186   }
 187
 188   TH1::AddDirectory(oldStatus);
 189
 190   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 191 }
 192
 193 //________________________________________________________________________
 194 Bool_t AliAnalysisTaskEmcalDiJetAna::FillHistograms()
 195 {
 196   // Fill histograms.
 197
 198   fh2CentRhoCh->Fill(fCent,fRhoChVal);
 199   fh2CentRhoScaled->Fill(fCent,fRhoFullVal);
 200
 201
 202   Int_t nJetsFull = GetNJets(fContainerFull);
 203   for(Int_t ij=0; ij<nJetsFull; ij++) {
 204     AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerFull));
 205     if(!jet) continue; //jet not selected
 206
 207     Double_t jetPt = GetJetPt(jet,0);
 208     fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 209
 210   }
 211
 212   Int_t nJetsCh = GetNJets(fContainerCharged);
 213   for(Int_t ij=0; ij<nJetsCh; ij++) {
 214     AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerCharged));
 215     if(!jet) continue; //jet not selected
 216
 217       Double_t jetPt = GetJetPt(jet,1);
 218       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 219
 220   }
 221
 222
 223   return kTRUE;
 224 }
 225
 226 //________________________________________________________________________
 227 Bool_t AliAnalysisTaskEmcalDiJetAna::Run()
 228 {
 229   // Run analysis code here, if needed. It will be executed before FillHistograms().
 230
 231   //Check if event is selected (vertex & pile-up)
 232   if(!SelectEvent())
 233     return kFALSE;
 234
 235   fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
 236
 237   if(fRhoType==0) {
 238     fRhoFullVal = 0.;
 239     fRhoChVal = 0.;
 240   }
 241   if(fRhoType==1) {
 242     fRhoFullVal = GetRhoVal(fContainerFull);
 243     fRhoChVal = GetRhoVal(fContainerCharged);
 244   }
 245
 246   // MatchFullAndChargedJets();
 247   if(fDoChargedCharged)   CorrelateJets(1);
 248
 249   if(fDoMatchFullCharged) FillMatchFullChargedHistos(fContainerFull,fContainerCharged);
 250
 251   if(fDoFullCharged)      {
 252     SetChargedFractionIndex();
 253     CorrelateJets(2);
 254   }
 255
 256   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 257 }
 258
 259 //________________________________________________________________________
 260 void AliAnalysisTaskEmcalDiJetAna::CorrelateJets(const Int_t type) {
 261   //
 262   // Correlate jets and fill histos
 263   //
 264
 265   Int_t typet = 0;
 266   Int_t typea = 0;
 267   if(type==0) { //full-full
 268     typet = fContainerFull;
 269     typea = fContainerFull;
 270   }
 271   else if(type==1) { //charged-charged
 272     typet = fContainerCharged;
 273     typea = fContainerCharged;
 274   }
 275   else if(type==2) { //full-charged
 276     typet = fContainerFull;
 277     typea = fContainerCharged;
 278   }
 279   else {
 280     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 281     return;
 282   }
 283
 284   Int_t nJetsTrig  = 0;
 285   Int_t nJetsAssoc = 0;
 286   if(type==0) {
 287     nJetsTrig  = GetNJets(fContainerFull);
 288     nJetsAssoc = nJetsTrig;
 289   }
 290   else if(type==1) {
 291     nJetsTrig  = GetNJets(fContainerCharged);
 292     nJetsAssoc = nJetsTrig;
 293   }
 294   else if(type==2) {
 295     nJetsTrig  = GetNJets(fContainerFull);
 296     nJetsAssoc = GetNJets(fContainerCharged);
 297   }
 298
 299   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 300
 301     AliEmcalJet *jetTrig = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typet));
 302     if(!jetTrig)
 303       continue; //jet not selected
 304
 305     Double_t jetTrigPt = GetJetPt(jetTrig,typet);
 306
 307     if(jetTrigPt<fPtMinTriggerJet)
 308       continue;
 309
 310     for(Int_t ija=0; ija<nJetsAssoc; ija++) {
 311       if(IsSameJet(ijt,ija,type)) continue;
 312
 313       AliEmcalJet *jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
 314       if(!jetAssoc)
 315         continue;
 316
 317       Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
 318
 319       if(jetTrigPt>jetAssocPt)
 320         FillDiJetHistos(jetTrig,jetAssoc, type);
 321
 322     } // associate jet loop
 323   }//trigger jet loop
 324
 325 }
 326
 327 //________________________________________________________________________
 328 void AliAnalysisTaskEmcalDiJetAna::FillDiJetHistos(const AliEmcalJet *jet1, const AliEmcalJet *jet2, const Int_t mode) {
 329   //
 330   // Fill histos
 331   // mode: charged+neutral = 0
 332   //       charged only    = 1
 333   //       full vs charged = 2
 334   //
 335
 336   Int_t typet = 0;
 337   Int_t typea = 0;
 338   if(mode==0) {
 339     typet = 0;
 340     typea = 0;
 341   }
 342   else if(mode==1) {
 343     typet = 1;
 344     typea = 1;
 345   }
 346   else if(mode==2) {
 347     typet = 0;
 348     typea = 1;
 349   }
 350   else {
 351     AliWarning(Form("%s: mode %d of dijet correlation not defined!",GetName(),mode));
 352     return;
 353   }
 354
 355   Double_t jetTrigPt = GetJetPt(jet1,typet);
 356   Double_t jetAssocPt = GetJetPt(jet2,typea);
 357
 358   Double_t deltaPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
 359   if(fDebug>10) Printf("deltaPhi:%.2f",deltaPhi);
 360
 361   Double_t kT = jetTrigPt*TMath::Sin(deltaPhi);
 362
 363   Double_t dijetEta = (jet1->Eta()+jet2->Eta())/2.;
 364
 365
 366   Double_t diJetVars[6] = {jetTrigPt,jetAssocPt,deltaPhi,kT,dijetEta,fCent};
 367
 368   if(mode==0)
 369     fhnDiJetVarsFull->Fill(diJetVars);
 370   else if(mode==1)
 371     fhnDiJetVarsCh->Fill(diJetVars);
 372   else if(mode==2)
 373     fhnDiJetVarsFullCharged->Fill(diJetVars);
 374
 375 }
 376
 377 //________________________________________________________________________
 378 void AliAnalysisTaskEmcalDiJetAna::FillMatchFullChargedHistos(Int_t cFull,Int_t cCharged) {
 379   //
 380   // Match full to charged jets and fill histo
 381   //
 382
 383   Int_t match = MatchFullAndChargedJets(cFull,cCharged);
 384   if(match==0) {
 385     if(fDebug>1) AliWarning(Form("%s: matching failed",GetName()));
 386     return;
 387   }
 388
 389   for(int ig = 0;ig < GetNJets(cFull);++ig){
 390     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ig, cFull));
 391     if(!jetFull) continue;
 392
 393     AliEmcalJet *jetCh = jetFull->ClosestJet();
 394     if(!jetCh) continue;
 395
 396     Double_t shFraction = GetFractionSharedPt(jetFull,jetCh);
 397     if(fDebug>10) Printf("shared charged pT:%.2f",shFraction);
 398     Double_t matchVars[7] = {
 399       jetFull->Pt(),
 400       jetCh->Pt(),
 401       GetDeltaPhi(jetFull->Phi(),jetCh->Phi()),
 402       jetFull->Eta()-jetCh->Eta(),
 403       GetDeltaR(jetFull,jetCh),
 404       shFraction,TMath::Min((Float_t)ig+0.5,2.5)
 405     };
 406       fhnMatchingFullCharged->Fill(matchVars);
 407
 408   }//loop over full jets
 409
 410 }
 411
 412
 413 //________________________________________________________________________
 414 Int_t AliAnalysisTaskEmcalDiJetAna::MatchFullAndChargedJets(Int_t cFull, Int_t cCharged) {
 415   //
 416   // Match charged jets to full jets
 417   //
 418
 419   if(GetNJets(cFull)<1) {
 420     if(fDebug>1) AliInfo(Form("%s: no full jets: %d", GetName(),GetNJets(cFull)));
 421     return 0;
 422   }
 423
 424   if(GetNJets(cCharged)<1) {
 425     if(fDebug>1) AliInfo(Form("%s: no charged jets: %d", GetName(),GetNJets(cCharged)));
 426     return 0;
 427   }
 428
 429   TClonesArray *cJetsFull = GetJetArray(cFull);
 430   TClonesArray *cJetsCharged = GetJetArray(cCharged);
 431
 432   if(!cJetsFull) {
 433     if(fDebug>1) AliInfo(Form("%s: no full jet array",GetName()));
 434     return 0;
 435   }
 436
 437   if(!cJetsCharged) {
 438     if(fDebug>1)  AliInfo(Form("%s: no charged jet array",GetName()));
 439     return 0;
 440   }
 441
 442
 443   if(!fMatchingDone) {
 444       MatchJetsGeo(cFull, cCharged, fDebug);
 445       return 1;
 446   } else {
 447     if(fDebug>1) AliInfo(Form("%s: Matching already done before",GetName()));
 448     return 1;
 449   }
 450
 451 }
 452
 453 //_______________________________________________________________________
 454 void AliAnalysisTaskEmcalDiJetAna::Terminate(Option_t *)
 455 {
 456   // Called once at the end of the analysis.
 457 }
 458
 459
 460