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   fh3AssocLostPtDeltaPhiCharged(0),
  46   fh3AssocLostPtDeltaPhiFull(0),
  47   fhnMatchingCharged(0),
  48   fhnMatchingFull(0),
  49   fnUsedResponseVar(0)
  50 {
  51   // Default constructor.
  52
  53   SetMakeGeneralHistograms(kTRUE);
  54
  55   for(Int_t i = 0; i<2; i++) {
  56     fh1TriggersCharged[i] = 0;
  57     fh1TriggersFull[i] = 0;
  58   }
  59
  60 }
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse(const char *name) :
  64   AliAnalysisTaskEmcalDiJetBase(name),
  65   fDoMatchFullCharged(kTRUE),
  66   fhnDiJetResponseCharged(0),
  67   fhnDiJetResponseFullCharged(0),
  68   fh1TriggersLostCharged(0),
  69   fh1TriggersLostFull(0),
  70   fh3AssocLostPtDeltaPhiCharged(0),
  71   fh3AssocLostPtDeltaPhiFull(0),
  72   fhnMatchingCharged(0),
  73   fhnMatchingFull(0),
  74   fnUsedResponseVar(0)
  75 {
  76   // Standard constructor.
  77
  78   SetMakeGeneralHistograms(kTRUE);
  79
  80   for(Int_t i = 0; i<2; i++) {
  81     fh1TriggersCharged[i] = 0;
  82     fh1TriggersFull[i] = 0;
  83   }
  84
  85 }
  86
  87 //________________________________________________________________________
  88 AliAnalysisTaskEmcalDiJetResponse::~AliAnalysisTaskEmcalDiJetResponse()
  89 {
  90   // Destructor.
  91 }
  92
  93
  94 //________________________________________________________________________
  95 void AliAnalysisTaskEmcalDiJetResponse::UserCreateOutputObjects()
  96 {
  97   // Create user output.
  98
  99   InitOnce();
 100
 101   AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
 102
 103   Bool_t oldStatus = TH1::AddDirectoryStatus();
 104   TH1::AddDirectory(kFALSE);
 105
 106   //Store dijet vars: pt,trig MC, pt,trig DET, pt,ass MC, pt,ass DET, dPhi MC, dPhi Det, kT MC, kT Det
 107   const Int_t nBinsSparse0 = 10;
 108   const Int_t nBinsPt = 250;
 109   const Int_t nBinsDPhi     = 36;
 110   const Int_t nBinsKt       = 25;
 111   const Int_t nBinsDiJetEta = 40;
 112   const Int_t nBinsAj       = 50;
 113   const Int_t nBinsVar[2] = {nBinsKt,nBinsDiJetEta};
 114
 115   const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsPt,nBinsPt,nBinsDPhi,nBinsDPhi,nBinsVar[fnUsedResponseVar],nBinsVar[fnUsedResponseVar],nBinsAj,nBinsAj};
 116
 117   const Double_t minPt = 0.;
 118   const Double_t maxPt = 250.;
 119   const Double_t minVar[2] = {   0.,-1.};
 120   const Double_t maxVar[2] = { 100., 1.};
 121
 122   const Double_t xmin0[nBinsSparse0]  = {  minPt, minPt, minPt, minPt, 0.5*TMath::Pi(), 0.5*TMath::Pi(), minVar[fnUsedResponseVar], minVar[fnUsedResponseVar],0.,0.};
 123   const Double_t xmax0[nBinsSparse0]  = {  maxPt, maxPt, maxPt, maxPt, 1.5*TMath::Pi(), 1.5*TMath::Pi(), maxVar[fnUsedResponseVar], maxVar[fnUsedResponseVar],1.,1.};
 124
 125   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};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
 126
 127   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};A_{j}^{part}A_{j}^{det}",nBinsSparse0,nBins0,xmin0,xmax0);
 128
 129   if(fnUsedResponseVar==1) {
 130     fhnDiJetResponseCharged->SetTitle("fhnDiJetResponseCharged DiJetEta");
 131     fhnDiJetResponseCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
 132     fhnDiJetResponseCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
 133
 134     fhnDiJetResponseFullCharged->SetTitle("fhnDiJetResponseFullCharged DiJetEta");
 135     fhnDiJetResponseFullCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
 136     fhnDiJetResponseFullCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
 137   }
 138
 139   fOutput->Add(fhnDiJetResponseCharged);
 140   fOutput->Add(fhnDiJetResponseFullCharged);
 141
 142   TString strType = "";
 143   for(Int_t i = 0; i<2; i++) {
 144     if(i==0)      strType="Part";
 145     else if(i==1) strType="Det";
 146     fh1TriggersCharged[i] = new TH1F(Form("fh1TriggersCharged%s",strType.Data()),Form("fh1TriggersCharged%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
 147     fOutput->Add(fh1TriggersCharged[i]);
 148
 149     fh1TriggersFull[i] = new TH1F(Form("fh1TriggersFull%s",strType.Data()),Form("fh1TriggersFull%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
 150     fOutput->Add(fh1TriggersFull[i]);
 151   }
 152
 153   fh1TriggersLostCharged = new TH1F("fh1TriggersLostCharged","fh1TriggersLostCharged;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
 154   fOutput->Add(fh1TriggersLostCharged);
 155
 156   fh1TriggersLostFull = new TH1F("fh1TriggersLostFull","fh1TriggersLostFull;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
 157   fOutput->Add(fh1TriggersLostFull);
 158
 159   fh3AssocLostPtDeltaPhiCharged = new TH3F("fh3AssocLostPtDeltaPhiCharged","fh3AssocLostPtDeltaPhiCharged;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
 160   fOutput->Add(fh3AssocLostPtDeltaPhiCharged);
 161
 162   fh3AssocLostPtDeltaPhiFull = new TH3F("fh3AssocLostPtDeltaPhiFull","fh3AssocLostPtDeltaPhiFull;p_{T,trig}^{ch};p_{T,assoc}^{ch};#Delta#varphi",nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
 163   fOutput->Add(fh3AssocLostPtDeltaPhiFull);
 164
 165   const Int_t nBinsSparseMatch = 6;
 166   const Int_t nBinsDPhiMatch = 80;
 167   const Int_t nBinsDEtaMatch = 80;
 168   const Int_t nBinsDR = 20;
 169   const Int_t nBinsType = 3;
 170   const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsType};
 171   //pTpart, pTdet, deltaPhi, deltaEta, deltaR, jet type (leading,subleading,other)
 172   const Double_t xminMatch[nBinsSparseMatch]  = { minPt, minPt, -0.5,-0.5, 0., 0};
 173   const Double_t xmaxMatch[nBinsSparseMatch]  = { maxPt, maxPt,  0.5, 0.5, 0.5,3};
 174   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",
 175                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 176   fOutput->Add(fhnMatchingCharged);
 177
 178   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",
 179                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 180   fOutput->Add(fhnMatchingFull);
 181
 182
 183   // =========== Switch on Sumw2 for all histos ===========
 184   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 185     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 186     if (h1){
 187       h1->Sumw2();
 188       continue;
 189     }
 190     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 191     if(hn)hn->Sumw2();
 192   }
 193
 194   TH1::AddDirectory(oldStatus);
 195
 196   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 197 }
 198
 199 //________________________________________________________________________
 200 Bool_t AliAnalysisTaskEmcalDiJetResponse::FillHistograms()
 201 {
 202   // Fill histograms.
 203
 204
 205
 206   return kTRUE;
 207 }
 208
 209 //________________________________________________________________________
 210 Bool_t AliAnalysisTaskEmcalDiJetResponse::Run()
 211 {
 212   // Run analysis code here, if needed. It will be executed before FillHistograms().
 213
 214   //Check if event is selected (vertex & pile-up)
 215   if(!SelectEvent())
 216     return kFALSE;
 217
 218   fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
 219
 220   if(fRhoType==0) {
 221     fRhoFullVal = 0.;
 222     fRhoChVal = 0.;
 223   }
 224   if(fRhoType==1) {
 225     fRhoFullVal = GetRhoVal(fContainerFull);
 226     fRhoChVal = GetRhoVal(fContainerCharged);
 227   }
 228
 229
 230   //Do matching
 231   MatchJetsGeo(fContainerCharged,fContainerChargedMC,0,0.3,1);
 232   MatchJetsGeo(fContainerFull,fContainerFullMC,0,0.3,2);
 233
 234   //Fill particle-detector level matching histos
 235
 236   if(fDoChargedCharged)   CorrelateJets(1);
 237
 238   if(fDoFullCharged) {
 239     SetChargedFractionIndexMC();
 240     CorrelateJets(2);
 241   }
 242
 243   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 244 }
 245
 246 //________________________________________________________________________
 247 void AliAnalysisTaskEmcalDiJetResponse::CorrelateJets(const Int_t type) {
 248   //
 249   // Correlate jets and fill histos
 250   //
 251
 252   if( fJetCorrelationType==kCorrelateAll )
 253     CorrelateAllJets(type);
 254   else if( fJetCorrelationType==kCorrelateTwo )
 255     CorrelateTwoJets(type);
 256   else if( fJetCorrelationType==kCorrelateLS )
 257     AliWarning(Form("%s: leading-subleading correlation not implemented for response!",GetName()));
 258
 259   return;
 260
 261 }
 262
 263
 264 //________________________________________________________________________
 265 void AliAnalysisTaskEmcalDiJetResponse::CorrelateAllJets(const Int_t type) {
 266   //
 267   // Correlate jets and fill histos
 268   //
 269
 270   Int_t typet = 0;
 271   Int_t typea = 0;
 272   Int_t typetMC = 0;
 273   Int_t typeaMC = 0;
 274   if(type==0) { //full-full
 275     typetMC = fContainerFullMC;
 276     typeaMC = fContainerFullMC;
 277     typet = fContainerFull;
 278     typea = fContainerFull;
 279   }
 280   else if(type==1) { //charged-charged
 281     typetMC = fContainerChargedMC;
 282     typeaMC = fContainerChargedMC;
 283     typet = fContainerCharged;
 284     typea = fContainerCharged;
 285   }
 286   else if(type==2) { //full-charged
 287     typetMC = fContainerFullMC;
 288     typeaMC = fContainerChargedMC;
 289     typet = fContainerFull;
 290     typea = fContainerCharged;
 291   }
 292   else {
 293     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 294     return;
 295   }
 296
 297   Int_t nJetsTrig  = 0;
 298   Int_t nJetsAssoc = 0;
 299   if(type==0) {
 300     nJetsTrig  = GetNJets(fContainerFullMC);
 301     nJetsAssoc = nJetsTrig;
 302   }
 303   else if(type==1) {
 304     nJetsTrig  = GetNJets(fContainerChargedMC);
 305     nJetsAssoc = nJetsTrig;
 306   }
 307   else if(type==2) {
 308     nJetsTrig  = GetNJets(fContainerFullMC);
 309     nJetsAssoc = GetNJets(fContainerChargedMC);
 310   }
 311
 312
 313   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 314
 315     AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
 316     if(!jetTrigMC) continue; //jet not selected
 317
 318     Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
 319
 320     if(jetTrigPtMC<fPtMinTriggerJet)
 321       continue;
 322
 323     if(type==1)
 324       fh1TriggersCharged[0]->Fill(jetTrigPtMC);
 325     if(type==2)
 326       fh1TriggersFull[0]->Fill(jetTrigPtMC);
 327
 328     AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
 329     if(!jetTrigDet) {
 330       //trigger is lost
 331       if(type==1)
 332         fh1TriggersLostCharged->Fill(jetTrigPtMC);
 333       if(type==2)
 334         fh1TriggersLostFull->Fill(jetTrigPtMC);
 335
 336       continue;
 337     }
 338
 339     if(type==1)
 340       fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
 341     if(type==2)
 342       fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
 343
 344     for(Int_t ija=0; ija<nJetsAssoc; ija++) {
 345       if(IsSameJet(ijt,ija,type,kTRUE)) continue;
 346
 347       AliEmcalJet *jetAssocMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typeaMC));
 348       if(!jetAssocMC) continue;
 349
 350       Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
 351
 352       //Now check if jets are also there on detector level
 353       AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
 354       if(!jetAssocDet) {
 355         //dijet is lost
 356       if(type==1)
 357         fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
 358       if(type==2)
 359         fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
 360         continue;
 361       }
 362
 363       FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
 364
 365     } // associate jet loop
 366   }//trigger jet loop
 367
 368 }
 369
 370 //________________________________________________________________________
 371 void AliAnalysisTaskEmcalDiJetResponse::CorrelateTwoJets(const Int_t type) {
 372   //
 373   // Correlate jets and fill histos
 374   //
 375
 376   Int_t typet = 0;
 377   Int_t typea = 0;
 378   Int_t typetMC = 0;
 379   Int_t typeaMC = 0;
 380   if(type==0) { //full-full
 381     typetMC = fContainerFullMC;
 382     typeaMC = fContainerFullMC;
 383     typet = fContainerFull;
 384     typea = fContainerFull;
 385   }
 386   else if(type==1) { //charged-charged
 387     typetMC = fContainerChargedMC;
 388     typeaMC = fContainerChargedMC;
 389     typet = fContainerCharged;
 390     typea = fContainerCharged;
 391   }
 392   else if(type==2) { //full-charged
 393     typetMC = fContainerFullMC;
 394     typeaMC = fContainerChargedMC;
 395     typet = fContainerFull;
 396     typea = fContainerCharged;
 397   }
 398   else {
 399     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 400     return;
 401   }
 402
 403   Int_t nJetsTrig  = 0;
 404   Int_t nJetsAssoc = 0;
 405   if(type==0) {
 406     nJetsTrig  = GetNJets(fContainerFullMC);
 407     nJetsAssoc = nJetsTrig;
 408   }
 409   else if(type==1) {
 410     nJetsTrig  = GetNJets(fContainerChargedMC);
 411     nJetsAssoc = nJetsTrig;
 412   }
 413   else if(type==2) {
 414     nJetsTrig  = GetNJets(fContainerFullMC);
 415     nJetsAssoc = GetNJets(fContainerChargedMC);
 416   }
 417
 418
 419   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 420
 421     AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
 422     if(!jetTrigMC) continue; //jet not selected
 423
 424     Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
 425
 426     if(jetTrigPtMC<fPtMinTriggerJet)
 427       continue;
 428
 429     if(type==1)
 430       fh1TriggersCharged[0]->Fill(jetTrigPtMC);
 431     if(type==2)
 432       fh1TriggersFull[0]->Fill(jetTrigPtMC);
 433
 434     AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
 435     if(!jetTrigDet) {
 436       //trigger is lost
 437       if(type==1)
 438         fh1TriggersLostCharged->Fill(jetTrigPtMC);
 439       if(type==2)
 440         fh1TriggersLostFull->Fill(jetTrigPtMC);
 441
 442       continue;
 443     }
 444
 445     if(type==1)
 446       fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
 447     if(type==2)
 448       fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
 449
 450
 451     AliEmcalJet *jetAssocMC = GetLeadingJetOppositeHemisphere(type,typeaMC,jetTrigMC);
 452     if(!jetAssocMC) continue;
 453
 454     Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
 455
 456     //Now check if jets are also there on detector level
 457     AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
 458     if(!jetAssocDet) {
 459       //dijet is lost
 460       if(type==1)
 461         fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
 462       if(type==2)
 463         fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
 464       continue;
 465     }
 466
 467     if(type==0 || type==1) {
 468       if(GetJetPt(jetAssocDet,typea)>GetJetPt(jetTrigDet,typet))
 469         continue;
 470     }
 471
 472     FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
 473
 474
 475   }//trigger jet loop
 476
 477 }
 478
 479 //________________________________________________________________________
 480 void AliAnalysisTaskEmcalDiJetResponse::FillDiJetResponse(const AliEmcalJet *jetTrigMC, const AliEmcalJet *jetAssocMC, const AliEmcalJet *jetTrigDet, const AliEmcalJet *jetAssocDet, Int_t type) {
 481
 482   //Fill dijet response
 483
 484   Int_t typet = 0;
 485   Int_t typea = 0;
 486   Int_t typetMC = 0;
 487   Int_t typeaMC = 0;
 488   if(type==0) { //full-full
 489     typetMC = fContainerFullMC;
 490     typeaMC = fContainerFullMC;
 491     typet = fContainerFull;
 492     typea = fContainerFull;
 493   }
 494   else if(type==1) { //charged-charged
 495     typetMC = fContainerChargedMC;
 496     typeaMC = fContainerChargedMC;
 497     typet = fContainerCharged;
 498     typea = fContainerCharged;
 499   }
 500   else if(type==2) { //full-charged
 501     typetMC = fContainerFullMC;
 502     typeaMC = fContainerChargedMC;
 503     typet = fContainerFull;
 504     typea = fContainerCharged;
 505   }
 506   else {
 507     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 508     return;
 509   }
 510
 511   Double_t jetTrigPtMC  = GetJetPt(jetTrigMC,typetMC);
 512   Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
 513
 514   Double_t varDet[2] = {TMath::Abs(GetJetPt(jetTrigDet,typet)*TMath::Sin(GetDeltaPhi(jetTrigDet,jetAssocDet))),(jetTrigDet->Eta()+jetAssocDet->Eta())/2.};
 515   Double_t varPart[2] = {TMath::Abs(jetTrigPtMC*TMath::Sin(GetDeltaPhi(jetTrigMC,jetAssocMC))),(jetTrigMC->Eta()+jetAssocMC->Eta())/2.};
 516
 517   Double_t ajDet  = (GetJetPt(jetTrigDet,typet)-GetJetPt(jetAssocDet,typea))/(GetJetPt(jetTrigDet,typet)+GetJetPt(jetAssocDet,typea));
 518   Double_t ajPart = (jetTrigPtMC-jetAssocPtMC)/(jetTrigPtMC+jetAssocPtMC);
 519
 520   //Store dijet vars: pt,trig MC; pt,trig DET; pt,ass MC; pt,ass DET; dPhi MC; dPhi Det; kT MC; kT Det;
 521   Double_t diJetVars[10] = {
 522     jetTrigPtMC,
 523     GetJetPt(jetTrigDet,typet),
 524     jetAssocPtMC,
 525     GetJetPt(jetAssocDet,typea),
 526     GetDeltaPhi(jetTrigMC,jetAssocMC),
 527     GetDeltaPhi(jetTrigDet,jetAssocDet),
 528     varPart[fnUsedResponseVar],
 529     varDet[fnUsedResponseVar],
 530     ajDet,
 531     ajPart
 532   };
 533
 534   if(type==1)
 535     fhnDiJetResponseCharged->Fill(diJetVars);
 536   else if(type==2)
 537     fhnDiJetResponseFullCharged->Fill(diJetVars);
 538
 539
 540 }
 541
 542 //________________________________________________________________________
 543 void AliAnalysisTaskEmcalDiJetResponse::FillMatchHistos() {
 544   //
 545   // Fill Particle-Detector level matching histos
 546   //
 547
 548   for(int i = 0; i < GetNJets(fContainerFull);++i) {
 549     AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerFull));
 550     if(!jetDet) continue;
 551
 552     AliEmcalJet *jetPart = jetDet->ClosestJet();
 553     if(!jetPart) continue;
 554
 555     Double_t matchVars[6] = {
 556       jetPart->Pt(),
 557       jetDet->Pt(),
 558       GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
 559       jetPart->Eta()-jetDet->Eta(),
 560       GetDeltaR(jetPart,jetDet),
 561       TMath::Min((Float_t)i+0.5,2.5)
 562     };
 563     fhnMatchingFull->Fill(matchVars);
 564
 565   }//loop over full jets
 566
 567   for(int i = 0; i < GetNJets(fContainerCharged);++i) {
 568     AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerCharged));
 569     if(!jetDet) continue;
 570
 571     AliEmcalJet *jetPart = jetDet->ClosestJet();
 572     if(!jetPart) continue;
 573
 574     Double_t matchVars[6] = {
 575       jetPart->Pt(),
 576       jetDet->Pt(),
 577       GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
 578       jetPart->Eta()-jetDet->Eta(),
 579       GetDeltaR(jetPart,jetDet),
 580       TMath::Min((Float_t)i+0.5,2.5)
 581     };
 582     fhnMatchingCharged->Fill(matchVars);
 583
 584   }//loop over charged jets
 585
 586 }
 587
 588
 589 //________________________________________________________________________
 590 Bool_t AliAnalysisTaskEmcalDiJetResponse::RetrieveEventObjects() {
 591   //
 592   // retrieve event objects
 593   //
 594
 595   if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
 596     return kFALSE;
 597
 598   return kTRUE;
 599
 600 }
 601
 602 //_______________________________________________________________________
 603 void AliAnalysisTaskEmcalDiJetResponse::Terminate(Option_t *)
 604 {
 605   // Called once at the end of the analysis.
 606 }