PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetResponse.cxx

   1 //
   2 // Dijet response 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 "AliAnalysisTaskEmcalDiJetResponse.h"
  34
  35 ClassImp(AliAnalysisTaskEmcalDiJetResponse)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse() :
  39   AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetResponse"),
  40   fDoMatchFullCharged(kTRUE),
  41   fhnDiJetResponseCharged(0),
  42   fhnDiJetResponseFullCharged(0),
  43   fh1TriggersLostCharged(0),
  44   fh1TriggersLostFull(0),
  45   fhnMatchingCharged(0),
  46   fhnMatchingFull(0)
  47 {
  48   // Default constructor.
  49
  50
  51   SetMakeGeneralHistograms(kTRUE);
  52
  53   for(Int_t i = 0; i<2; i++) {
  54     fh1TriggersCharged[i] = 0;
  55     fh1TriggersFull[i] = 0;
  56   }
  57
  58 }
  59
  60 //________________________________________________________________________
  61 AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse(const char *name) :
  62   AliAnalysisTaskEmcalDiJetBase(name),
  63   fDoMatchFullCharged(kTRUE),
  64   fhnDiJetResponseCharged(0),
  65   fhnDiJetResponseFullCharged(0),
  66   fh1TriggersLostCharged(0),
  67   fh1TriggersLostFull(0),
  68   fhnMatchingCharged(0),
  69   fhnMatchingFull(0)
  70 {
  71   // Standard constructor.
  72
  73   SetMakeGeneralHistograms(kTRUE);
  74
  75   for(Int_t i = 0; i<2; i++) {
  76     fh1TriggersCharged[i] = 0;
  77     fh1TriggersFull[i] = 0;
  78   }
  79
  80 }
  81
  82 //________________________________________________________________________
  83 AliAnalysisTaskEmcalDiJetResponse::~AliAnalysisTaskEmcalDiJetResponse()
  84 {
  85   // Destructor.
  86 }
  87
  88
  89 //________________________________________________________________________
  90 void AliAnalysisTaskEmcalDiJetResponse::UserCreateOutputObjects()
  91 {
  92   // Create user output.
  93
  94   InitOnce();
  95
  96   AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
  97
  98   Bool_t oldStatus = TH1::AddDirectoryStatus();
  99   TH1::AddDirectory(kFALSE);
 100
 101   //Store dijet vars: pt,trig MC, pt,trig DET, pt,ass MC, pt,ass DET, dPhi MC, dPhi Det, kT MC, kT Det
 102   const Int_t nBinsSparse0 = 8;
 103   const Int_t nBinsPt = 250;
 104   const Int_t nBinsDPhi     = 72;
 105   const Int_t nBinsKt       = 50;
 106   const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsPt,nBinsPt,nBinsDPhi,nBinsDPhi,nBinsKt,nBinsKt};
 107
 108   const Double_t minPt = 0.;
 109   const Double_t maxPt = 250.;
 110
 111   const Double_t xmin0[nBinsSparse0]  = {  minPt, minPt, minPt, minPt, -0.5*TMath::Pi(), -0.5*TMath::Pi(), -100.,-100.};
 112   const Double_t xmax0[nBinsSparse0]  = {  maxPt, maxPt, maxPt, maxPt,  1.5*TMath::Pi(), 1.5*TMath::Pi(),  100., 100.};
 113
 114   fhnDiJetResponseCharged = new THnSparseF("fhnDiJetResponseCharged","fhnDiJetResponseCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
 115   fOutput->Add(fhnDiJetResponseCharged);
 116
 117   fhnDiJetResponseFullCharged = new THnSparseF("fhnDiJetResponseFullCharged","fhnDiJetResponseFullCharged;p_{T,trig}^{part};p_{T,trig}^{det};p_{T,ass}^{part};p_{T,ass}^{det};#Delta#varphi_{part};#Delta#varphi_{det};k_{T}^{part},k_{T}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
 118   fOutput->Add(fhnDiJetResponseFullCharged);
 119
 120   TString strType = "";
 121   for(Int_t i = 0; i<2; i++) {
 122     if(i==0)      strType="Part";
 123     else if(i==1) strType="Det";
 124     fh1TriggersCharged[i] = new TH1F(Form("fh1TriggersCharged%s",strType.Data()),Form("fh1TriggersCharged%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
 125     fOutput->Add(fh1TriggersCharged[i]);
 126
 127     fh1TriggersFull[i] = new TH1F(Form("fh1TriggersFull%s",strType.Data()),Form("fh1TriggersFull%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
 128     fOutput->Add(fh1TriggersFull[i]);
 129   }
 130
 131   fh1TriggersLostCharged = new TH1F("fh1TriggersLostCharged","fh1TriggersLostCharged;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
 132   fOutput->Add(fh1TriggersLostCharged);
 133
 134   fh1TriggersLostFull = new TH1F("fh1TriggersLostFull","fh1TriggersLostFull;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
 135   fOutput->Add(fh1TriggersLostFull);
 136
 137   const Int_t nBinsSparseMatch = 6;
 138   const Int_t nBinsDPhiMatch = 80;
 139   const Int_t nBinsDEtaMatch = 80;
 140   const Int_t nBinsDR = 20;
 141   const Int_t nBinsType = 3;
 142   const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsType};
 143   //pTpart, pTdet, deltaPhi, deltaEta, deltaR, jet type (leading,subleading,other)
 144   const Double_t xminMatch[nBinsSparseMatch]  = { minPt, minPt, -0.5,-0.5, 0., 0};
 145   const Double_t xmaxMatch[nBinsSparseMatch]  = { maxPt, maxPt,  0.5, 0.5, 0.5,3};
 146   fhnMatchingCharged = new THnSparseF("fhnMatchingCharged","fhnMatchingCharged;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
 147                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 148   fOutput->Add(fhnMatchingCharged);
 149
 150   fhnMatchingFull = new THnSparseF("fhnMatchingFull","fhnMatchingFull;#it{p}_{T,part} (GeV/#it{c});#it{p}_{T,det} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;type",
 151                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 152   fOutput->Add(fhnMatchingFull);
 153
 154
 155   // =========== Switch on Sumw2 for all histos ===========
 156   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 157     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 158     if (h1){
 159       h1->Sumw2();
 160       continue;
 161     }
 162     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 163     if(hn)hn->Sumw2();
 164   }
 165
 166   TH1::AddDirectory(oldStatus);
 167
 168   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 169 }
 170
 171 //________________________________________________________________________
 172 Bool_t AliAnalysisTaskEmcalDiJetResponse::FillHistograms()
 173 {
 174   // Fill histograms.
 175
 176
 177
 178   return kTRUE;
 179 }
 180
 181 //________________________________________________________________________
 182 Bool_t AliAnalysisTaskEmcalDiJetResponse::Run()
 183 {
 184   // Run analysis code here, if needed. It will be executed before FillHistograms().
 185
 186   //Check if event is selected (vertex & pile-up)
 187   if(!SelectEvent())
 188     return kFALSE;
 189
 190   fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
 191
 192   if(fRhoType==0) {
 193     fRhoFullVal = 0.;
 194     fRhoChVal = 0.;
 195   }
 196   if(fRhoType==1) {
 197     fRhoFullVal = GetRhoVal(fContainerFull);
 198     fRhoChVal = GetRhoVal(fContainerCharged);
 199   }
 200
 201
 202   //Do matching
 203   MatchJetsGeo(fContainerCharged,fContainerChargedMC,0,0.3,1);
 204   MatchJetsGeo(fContainerFull,fContainerFullMC,0,0.3,2);
 205
 206   //Fill particle-detector level matching histos
 207
 208   if(fDoChargedCharged)   CorrelateJets(1);
 209
 210   if(fDoFullCharged) {
 211     SetChargedFractionIndexMC();
 212     CorrelateJets(2);
 213   }
 214
 215   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 216 }
 217
 218 //________________________________________________________________________
 219 void AliAnalysisTaskEmcalDiJetResponse::CorrelateJets(const Int_t type) {
 220   //
 221   // Correlate jets and fill histos
 222   //
 223
 224   Int_t typet = 0;
 225   Int_t typea = 0;
 226   Int_t typetMC = 0;
 227   Int_t typeaMC = 0;
 228   if(type==0) { //full-full
 229     typetMC = fContainerFullMC;
 230     typeaMC = fContainerFullMC;
 231     typet = fContainerFull;
 232     typea = fContainerFull;
 233   }
 234   else if(type==1) { //charged-charged
 235     typetMC = fContainerChargedMC;
 236     typeaMC = fContainerChargedMC;
 237     typet = fContainerCharged;
 238     typea = fContainerCharged;
 239   }
 240   else if(type==2) { //full-charged
 241     typetMC = fContainerFullMC;
 242     typeaMC = fContainerChargedMC;
 243     typet = fContainerFull;
 244     typea = fContainerCharged;
 245   }
 246   else {
 247     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 248     return;
 249   }
 250
 251   Int_t nJetsTrig  = 0;
 252   Int_t nJetsAssoc = 0;
 253   if(type==0) {
 254     nJetsTrig  = GetNJets(fContainerFullMC);
 255     nJetsAssoc = nJetsTrig;
 256   }
 257   else if(type==1) {
 258     nJetsTrig  = GetNJets(fContainerChargedMC);
 259     nJetsAssoc = nJetsTrig;
 260   }
 261   else if(type==2) {
 262     nJetsTrig  = GetNJets(fContainerFullMC);
 263     nJetsAssoc = GetNJets(fContainerChargedMC);
 264   }
 265
 266
 267   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 268
 269     AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
 270     if(!jetTrigMC) continue; //jet not selected
 271
 272     Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
 273
 274     if(jetTrigPtMC<fPtMinTriggerJet)
 275       continue;
 276
 277     if(type==1)
 278       fh1TriggersCharged[0]->Fill(jetTrigPtMC);
 279     if(type==2)
 280       fh1TriggersFull[0]->Fill(jetTrigPtMC);
 281
 282     AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
 283     if(!jetTrigDet) {
 284       //trigger is lost
 285       if(type==1)
 286         fh1TriggersLostCharged->Fill(jetTrigPtMC);
 287       if(type==2)
 288         fh1TriggersLostFull->Fill(jetTrigPtMC);
 289
 290       continue;
 291     }
 292
 293     if(type==1)
 294       fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
 295     if(type==2)
 296       fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
 297
 298     for(Int_t ija=0; ija<nJetsAssoc; ija++) {
 299       if(IsSameJet(ijt,ija,type,kTRUE)) continue;
 300
 301       AliEmcalJet *jetAssocMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typeaMC));
 302       if(!jetAssocMC) continue;
 303
 304       Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
 305
 306       //Now check if jets are also there on detector level
 307       AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
 308       if(!jetTrigDet || !jetAssocDet) {
 309         //dijet is lost
 310         continue;
 311       }
 312
 313       //Store dijet vars: pt,trig MC; pt,trig DET; pt,ass MC; pt,ass DET; dPhi MC; dPhi Det; kT MC; kT Det;
 314       Double_t diJetVars[8] = {
 315         jetTrigPtMC,
 316         GetJetPt(jetTrigDet,typet),
 317         jetAssocPtMC,
 318         GetJetPt(jetAssocDet,typea),
 319         GetDeltaPhi(jetTrigMC,jetAssocMC),
 320         GetDeltaPhi(jetTrigDet,jetAssocDet),
 321         jetTrigPtMC*TMath::Sin(GetDeltaPhi(jetTrigMC,jetAssocMC)),
 322         GetJetPt(jetTrigDet,typet)*TMath::Sin(GetDeltaPhi(jetTrigDet,jetAssocDet))
 323       };
 324
 325       if(type==1)
 326         fhnDiJetResponseCharged->Fill(diJetVars);
 327       else if(type==2)
 328         fhnDiJetResponseFullCharged->Fill(diJetVars);
 329
 330
 331     } // associate jet loop
 332   }//trigger jet loop
 333
 334 }
 335
 336 //________________________________________________________________________
 337 void AliAnalysisTaskEmcalDiJetResponse::FillMatchHistos() {
 338   //
 339   // Fill Particle-Detector level matching histos
 340   //
 341
 342   for(int i = 0; i < GetNJets(fContainerFull);++i) {
 343     AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerFull));
 344     if(!jetDet) continue;
 345
 346     AliEmcalJet *jetPart = jetDet->ClosestJet();
 347     if(!jetPart) continue;
 348
 349     Double_t matchVars[7] = {
 350       jetPart->Pt(),
 351       jetDet->Pt(),
 352       GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
 353       jetPart->Eta()-jetDet->Eta(),
 354       GetDeltaR(jetPart,jetDet),
 355       TMath::Min((Float_t)i+0.5,2.5)
 356     };
 357     fhnMatchingFull->Fill(matchVars);
 358
 359   }//loop over full jets
 360
 361   for(int i = 0; i < GetNJets(fContainerCharged);++i) {
 362     AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerCharged));
 363     if(!jetDet) continue;
 364
 365     AliEmcalJet *jetPart = jetDet->ClosestJet();
 366     if(!jetPart) continue;
 367
 368     Double_t matchVars[7] = {
 369       jetPart->Pt(),
 370       jetDet->Pt(),
 371       GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
 372       jetPart->Eta()-jetDet->Eta(),
 373       GetDeltaR(jetPart,jetDet),
 374       TMath::Min((Float_t)i+0.5,2.5)
 375     };
 376     fhnMatchingCharged->Fill(matchVars);
 377
 378   }//loop over charged jets
 379
 380 }
 381
 382
 383 //________________________________________________________________________
 384 Bool_t AliAnalysisTaskEmcalDiJetResponse::RetrieveEventObjects() {
 385   //
 386   // retrieve event objects
 387   //
 388
 389   if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
 390     return kFALSE;
 391
 392   return kTRUE;
 393
 394 }
 395
 396 //_______________________________________________________________________
 397 void AliAnalysisTaskEmcalDiJetResponse::Terminate(Option_t *)
 398 {
 399   // Called once at the end of the analysis.
 400 }