2 // Dijet response analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.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"
32 #include "AliAnalysisTaskEmcalDiJetResponse.h"
34 ClassImp(AliAnalysisTaskEmcalDiJetResponse)
36 //________________________________________________________________________
37 AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse() :
38 AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetResponse"),
39 fDoMatchFullCharged(kTRUE),
40 fhnDiJetResponseCharged(0),
41 fhnDiJetResponseFullCharged(0),
42 fh1TriggersLostCharged(0),
43 fh1TriggersLostFull(0),
44 fh3AssocLostPtDeltaPhiCharged(0),
45 fh3AssocLostPtDeltaPhiFull(0),
46 fhnMatchingCharged(0),
50 // Default constructor.
52 SetMakeGeneralHistograms(kTRUE);
54 for(Int_t i = 0; i<2; i++) {
55 fh1TriggersCharged[i] = 0;
56 fh1TriggersFull[i] = 0;
61 //________________________________________________________________________
62 AliAnalysisTaskEmcalDiJetResponse::AliAnalysisTaskEmcalDiJetResponse(const char *name) :
63 AliAnalysisTaskEmcalDiJetBase(name),
64 fDoMatchFullCharged(kTRUE),
65 fhnDiJetResponseCharged(0),
66 fhnDiJetResponseFullCharged(0),
67 fh1TriggersLostCharged(0),
68 fh1TriggersLostFull(0),
69 fh3AssocLostPtDeltaPhiCharged(0),
70 fh3AssocLostPtDeltaPhiFull(0),
71 fhnMatchingCharged(0),
75 // Standard constructor.
77 SetMakeGeneralHistograms(kTRUE);
79 for(Int_t i = 0; i<2; i++) {
80 fh1TriggersCharged[i] = 0;
81 fh1TriggersFull[i] = 0;
86 //________________________________________________________________________
87 AliAnalysisTaskEmcalDiJetResponse::~AliAnalysisTaskEmcalDiJetResponse()
93 //________________________________________________________________________
94 void AliAnalysisTaskEmcalDiJetResponse::UserCreateOutputObjects()
96 // Create user output.
98 AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
100 Bool_t oldStatus = TH1::AddDirectoryStatus();
101 TH1::AddDirectory(kFALSE);
103 //Store dijet vars: pt,trig MC, pt,trig DET, pt,ass MC, pt,ass DET, dPhi MC, dPhi Det, kT MC, kT Det
104 const Int_t nBinsSparse0 = 10;
105 const Int_t nBinsPt = 250;
106 const Int_t nBinsDPhi = 36;
107 const Int_t nBinsKt = 25;
108 const Int_t nBinsDiJetEta = 40;
109 const Int_t nBinsAj = 50;
110 const Int_t nBinsVar[2] = {nBinsKt,nBinsDiJetEta};
112 const Int_t nBins0[nBinsSparse0] = {nBinsPt,nBinsPt,nBinsPt,nBinsPt,nBinsDPhi,nBinsDPhi,nBinsVar[fnUsedResponseVar],nBinsVar[fnUsedResponseVar],nBinsAj,nBinsAj};
114 const Double_t minPt = 0.;
115 const Double_t maxPt = 250.;
116 const Double_t minVar[2] = { 0.,-1.};
117 const Double_t maxVar[2] = { 100., 1.};
119 const Double_t xmin0[nBinsSparse0] = { minPt, minPt, minPt, minPt, 0.5*TMath::Pi(), 0.5*TMath::Pi(), minVar[fnUsedResponseVar], minVar[fnUsedResponseVar],0.,0.};
120 const Double_t xmax0[nBinsSparse0] = { maxPt, maxPt, maxPt, maxPt, 1.5*TMath::Pi(), 1.5*TMath::Pi(), maxVar[fnUsedResponseVar], maxVar[fnUsedResponseVar],1.,1.};
122 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);
124 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);
126 if(fnUsedResponseVar==1) {
127 fhnDiJetResponseCharged->SetTitle("fhnDiJetResponseCharged DiJetEta");
128 fhnDiJetResponseCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
129 fhnDiJetResponseCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
131 fhnDiJetResponseFullCharged->SetTitle("fhnDiJetResponseFullCharged DiJetEta");
132 fhnDiJetResponseFullCharged->GetAxis(6)->SetTitle("#eta_{dijet}^{part}");
133 fhnDiJetResponseFullCharged->GetAxis(7)->SetTitle("#eta_{dijet}^{det}");
136 fOutput->Add(fhnDiJetResponseCharged);
137 fOutput->Add(fhnDiJetResponseFullCharged);
139 TString strType = "";
140 for(Int_t i = 0; i<2; i++) {
141 if(i==0) strType="Part";
142 else if(i==1) strType="Det";
143 fh1TriggersCharged[i] = new TH1F(Form("fh1TriggersCharged%s",strType.Data()),Form("fh1TriggersCharged%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
144 fOutput->Add(fh1TriggersCharged[i]);
146 fh1TriggersFull[i] = new TH1F(Form("fh1TriggersFull%s",strType.Data()),Form("fh1TriggersFull%s;p_{T,trig}^{ch}",strType.Data()),nBinsPt,minPt,maxPt);
147 fOutput->Add(fh1TriggersFull[i]);
150 fh1TriggersLostCharged = new TH1F("fh1TriggersLostCharged","fh1TriggersLostCharged;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
151 fOutput->Add(fh1TriggersLostCharged);
153 fh1TriggersLostFull = new TH1F("fh1TriggersLostFull","fh1TriggersLostFull;p_{T,trig}^{ch}",nBinsPt,minPt,maxPt);
154 fOutput->Add(fh1TriggersLostFull);
156 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());
157 fOutput->Add(fh3AssocLostPtDeltaPhiCharged);
159 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());
160 fOutput->Add(fh3AssocLostPtDeltaPhiFull);
162 const Int_t nBinsSparseMatch = 6;
163 const Int_t nBinsDPhiMatch = 80;
164 const Int_t nBinsDEtaMatch = 80;
165 const Int_t nBinsDR = 20;
166 const Int_t nBinsType = 3;
167 const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsType};
168 //pTpart, pTdet, deltaPhi, deltaEta, deltaR, jet type (leading,subleading,other)
169 const Double_t xminMatch[nBinsSparseMatch] = { minPt, minPt, -0.5,-0.5, 0., 0};
170 const Double_t xmaxMatch[nBinsSparseMatch] = { maxPt, maxPt, 0.5, 0.5, 0.5,3};
171 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",
172 nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
173 fOutput->Add(fhnMatchingCharged);
175 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",
176 nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
177 fOutput->Add(fhnMatchingFull);
180 // =========== Switch on Sumw2 for all histos ===========
181 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
182 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
187 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
191 TH1::AddDirectory(oldStatus);
193 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
196 //________________________________________________________________________
197 Bool_t AliAnalysisTaskEmcalDiJetResponse::Run()
199 // Run analysis code here, if needed. It will be executed before FillHistograms().
201 //Check if event is selected (vertex & pile-up)
210 fRhoFullVal = GetRhoVal(fContainerFull);
211 fRhoChVal = GetRhoVal(fContainerCharged);
216 MatchJetsGeo(fContainerCharged,fContainerChargedMC,0,0.3,1);
217 MatchJetsGeo(fContainerFull,fContainerFullMC,0,0.3,2);
219 //Fill particle-detector level matching histos
221 if(fDoChargedCharged) CorrelateJets(1);
224 SetChargedFractionIndexMC();
228 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
231 //________________________________________________________________________
232 void AliAnalysisTaskEmcalDiJetResponse::CorrelateJets(const Int_t type) {
234 // Correlate jets and fill histos
237 if( fJetCorrelationType==kCorrelateAll )
238 CorrelateAllJets(type);
239 else if( fJetCorrelationType==kCorrelateTwo )
240 CorrelateTwoJets(type);
241 else if( fJetCorrelationType==kCorrelateLS )
242 AliWarning(Form("%s: leading-subleading correlation not implemented for response!",GetName()));
249 //________________________________________________________________________
250 void AliAnalysisTaskEmcalDiJetResponse::CorrelateAllJets(const Int_t type) {
252 // Correlate jets and fill histos
258 if(type==0) { //full-full
259 typetMC = fContainerFullMC;
260 typeaMC = fContainerFullMC;
261 typet = fContainerFull;
263 else if(type==1) { //charged-charged
264 typetMC = fContainerChargedMC;
265 typeaMC = fContainerChargedMC;
266 typet = fContainerCharged;
268 else if(type==2) { //full-charged
269 typetMC = fContainerFullMC;
270 typeaMC = fContainerChargedMC;
271 typet = fContainerFull;
274 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
279 Int_t nJetsAssoc = 0;
281 nJetsTrig = GetNJets(fContainerFullMC);
282 nJetsAssoc = nJetsTrig;
285 nJetsTrig = GetNJets(fContainerChargedMC);
286 nJetsAssoc = nJetsTrig;
289 nJetsTrig = GetNJets(fContainerFullMC);
290 nJetsAssoc = GetNJets(fContainerChargedMC);
294 for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
296 AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
297 if(!jetTrigMC) continue; //jet not selected
299 Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
301 if(jetTrigPtMC<fPtMinTriggerJet)
305 fh1TriggersCharged[0]->Fill(jetTrigPtMC);
307 fh1TriggersFull[0]->Fill(jetTrigPtMC);
309 AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
313 fh1TriggersLostCharged->Fill(jetTrigPtMC);
315 fh1TriggersLostFull->Fill(jetTrigPtMC);
321 fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
323 fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
325 for(Int_t ija=0; ija<nJetsAssoc; ija++) {
326 if(IsSameJet(ijt,ija,type,kTRUE)) continue;
328 AliEmcalJet *jetAssocMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typeaMC));
329 if(!jetAssocMC) continue;
331 Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
333 //Now check if jets are also there on detector level
334 AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
338 fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
340 fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
344 FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
346 } // associate jet loop
351 //________________________________________________________________________
352 void AliAnalysisTaskEmcalDiJetResponse::CorrelateTwoJets(const Int_t type) {
354 // Correlate jets and fill histos
361 if(type==0) { //full-full
362 typetMC = fContainerFullMC;
363 typeaMC = fContainerFullMC;
364 typet = fContainerFull;
365 typea = fContainerFull;
367 else if(type==1) { //charged-charged
368 typetMC = fContainerChargedMC;
369 typeaMC = fContainerChargedMC;
370 typet = fContainerCharged;
371 typea = fContainerCharged;
373 else if(type==2) { //full-charged
374 typetMC = fContainerFullMC;
375 typeaMC = fContainerChargedMC;
376 typet = fContainerFull;
377 typea = fContainerCharged;
380 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
386 nJetsTrig = GetNJets(fContainerFullMC);
389 nJetsTrig = GetNJets(fContainerChargedMC);
392 nJetsTrig = GetNJets(fContainerFullMC);
395 for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
397 AliEmcalJet *jetTrigMC = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typetMC));
398 if(!jetTrigMC) continue; //jet not selected
400 Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
402 if(jetTrigPtMC<fPtMinTriggerJet)
406 fh1TriggersCharged[0]->Fill(jetTrigPtMC);
408 fh1TriggersFull[0]->Fill(jetTrigPtMC);
410 AliEmcalJet *jetTrigDet = jetTrigMC->ClosestJet();
414 fh1TriggersLostCharged->Fill(jetTrigPtMC);
416 fh1TriggersLostFull->Fill(jetTrigPtMC);
422 fh1TriggersCharged[1]->Fill(GetJetPt(jetTrigDet,typet));
424 fh1TriggersFull[1]->Fill(GetJetPt(jetTrigDet,typet));
427 AliEmcalJet *jetAssocMC = GetLeadingJetOppositeHemisphere(type,typeaMC,jetTrigMC);
428 if(!jetAssocMC) continue;
430 Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
432 //Now check if jets are also there on detector level
433 AliEmcalJet *jetAssocDet = jetAssocMC->ClosestJet();
437 fh3AssocLostPtDeltaPhiCharged->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
439 fh3AssocLostPtDeltaPhiFull->Fill(jetTrigPtMC,jetAssocPtMC,GetDeltaPhi(jetTrigMC,jetAssocMC));
443 if(type==0 || type==1) {
444 if(GetJetPt(jetAssocDet,typea)>GetJetPt(jetTrigDet,typet))
448 FillDiJetResponse(jetTrigMC,jetAssocMC,jetTrigDet,jetAssocDet,type);
455 //________________________________________________________________________
456 void AliAnalysisTaskEmcalDiJetResponse::FillDiJetResponse(const AliEmcalJet *jetTrigMC, const AliEmcalJet *jetAssocMC, const AliEmcalJet *jetTrigDet, const AliEmcalJet *jetAssocDet, Int_t type) {
458 //Fill dijet response
464 if(type==0) { //full-full
465 typetMC = fContainerFullMC;
466 typeaMC = fContainerFullMC;
467 typet = fContainerFull;
468 typea = fContainerFull;
470 else if(type==1) { //charged-charged
471 typetMC = fContainerChargedMC;
472 typeaMC = fContainerChargedMC;
473 typet = fContainerCharged;
474 typea = fContainerCharged;
476 else if(type==2) { //full-charged
477 typetMC = fContainerFullMC;
478 typeaMC = fContainerChargedMC;
479 typet = fContainerFull;
480 typea = fContainerCharged;
483 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
487 Double_t jetTrigPtMC = GetJetPt(jetTrigMC,typetMC);
488 Double_t jetAssocPtMC = GetJetPt(jetAssocMC,typeaMC);
490 Double_t varDet[2] = {TMath::Abs(GetJetPt(jetTrigDet,typet)*TMath::Sin(GetDeltaPhi(jetTrigDet,jetAssocDet))),(jetTrigDet->Eta()+jetAssocDet->Eta())/2.};
491 Double_t varPart[2] = {TMath::Abs(jetTrigPtMC*TMath::Sin(GetDeltaPhi(jetTrigMC,jetAssocMC))),(jetTrigMC->Eta()+jetAssocMC->Eta())/2.};
493 Double_t ajDet = (GetJetPt(jetTrigDet,typet)-GetJetPt(jetAssocDet,typea))/(GetJetPt(jetTrigDet,typet)+GetJetPt(jetAssocDet,typea));
494 Double_t ajPart = (jetTrigPtMC-jetAssocPtMC)/(jetTrigPtMC+jetAssocPtMC);
496 //Store dijet vars: pt,trig MC; pt,trig DET; pt,ass MC; pt,ass DET; dPhi MC; dPhi Det; kT MC; kT Det;
497 Double_t diJetVars[10] = {
499 GetJetPt(jetTrigDet,typet),
501 GetJetPt(jetAssocDet,typea),
502 GetDeltaPhi(jetTrigMC,jetAssocMC),
503 GetDeltaPhi(jetTrigDet,jetAssocDet),
504 varPart[fnUsedResponseVar],
505 varDet[fnUsedResponseVar],
511 fhnDiJetResponseCharged->Fill(diJetVars);
513 fhnDiJetResponseFullCharged->Fill(diJetVars);
518 //________________________________________________________________________
519 void AliAnalysisTaskEmcalDiJetResponse::FillMatchHistos() {
521 // Fill Particle-Detector level matching histos
524 for(int i = 0; i < GetNJets(fContainerFull);++i) {
525 AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerFull));
526 if(!jetDet) continue;
528 AliEmcalJet *jetPart = jetDet->ClosestJet();
529 if(!jetPart) continue;
531 Double_t matchVars[6] = {
534 GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
535 jetPart->Eta()-jetDet->Eta(),
536 GetDeltaR(jetPart,jetDet),
537 TMath::Min((Float_t)i+0.5,2.5)
539 fhnMatchingFull->Fill(matchVars);
541 }//loop over full jets
543 for(int i = 0; i < GetNJets(fContainerCharged);++i) {
544 AliEmcalJet *jetDet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerCharged));
545 if(!jetDet) continue;
547 AliEmcalJet *jetPart = jetDet->ClosestJet();
548 if(!jetPart) continue;
550 Double_t matchVars[6] = {
553 GetDeltaPhi(jetPart->Phi(),jetDet->Phi()),
554 jetPart->Eta()-jetDet->Eta(),
555 GetDeltaR(jetPart,jetDet),
556 TMath::Min((Float_t)i+0.5,2.5)
558 fhnMatchingCharged->Fill(matchVars);
560 }//loop over charged jets
565 //________________________________________________________________________
566 Bool_t AliAnalysisTaskEmcalDiJetResponse::RetrieveEventObjects() {
568 // retrieve event objects
571 if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
578 //_______________________________________________________________________
579 void AliAnalysisTaskEmcalDiJetResponse::Terminate(Option_t *)
581 // Called once at the end of the analysis.