2 // Jet mass background analysis task.
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
12 #include <TLorentzVector.h>
20 #include "AliVCluster.h"
21 #include "AliVTrack.h"
22 #include "AliEmcalJet.h"
23 #include "AliRhoParameter.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 #include "AliClusterContainer.h"
33 #include "AliParticleContainer.h"
35 #include "AliAnalysisTaskEmcalJetMassBkg.h"
37 ClassImp(AliAnalysisTaskEmcalJetMassBkg)
39 //________________________________________________________________________
40 AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg() :
41 AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassBkg", kTRUE),
49 fConeMaxPhi(TMath::Pi()*2),
55 fh2PtVsMassRCExLJDPhi(0),
57 fh2PtVsMassPerpConeLJ(0),
58 fpPtVsMassPerpConeLJ(0),
59 fh2PtVsMassPerpConeTJ(0),
60 fpPtVsMassPerpConeTJ(0),
62 fh2CentVsMassRCExLJ(0),
63 fh2CentVsMassPerpConeLJ(0),
64 fh2CentVsMassPerpConeTJ(0),
66 fh2MultVsMassRCExLJ(0),
67 fh2MultVsMassPerpConeLJ(0),
68 fh2MultVsMassPerpConeTJ(0)
70 // Default constructor.
72 fh2PtVsMassRC = new TH2F*[fNcentBins];
73 fpPtVsMassRC = new TProfile*[fNcentBins];
74 fh2PtVsMassRCExLJDPhi = new TH3F*[fNcentBins];
75 fpPtVsMassRCExLJ = new TProfile*[fNcentBins];
76 fh2PtVsMassPerpConeLJ = new TH2F*[fNcentBins];
77 fpPtVsMassPerpConeLJ = new TProfile*[fNcentBins];
78 fh2PtVsMassPerpConeTJ = new TH2F*[fNcentBins];
79 fpPtVsMassPerpConeTJ = new TProfile*[fNcentBins];
81 for (Int_t i = 0; i < fNcentBins; i++) {
84 fh2PtVsMassRCExLJDPhi[i] = 0;
85 fpPtVsMassRCExLJ[i] = 0;
86 fh2PtVsMassPerpConeLJ[i] = 0;
87 fpPtVsMassPerpConeLJ[i] = 0;
88 fh2PtVsMassPerpConeTJ[i] = 0;
89 fpPtVsMassPerpConeTJ[i] = 0;
92 SetMakeGeneralHistograms(kTRUE);
96 //________________________________________________________________________
97 AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg(const char *name) :
98 AliAnalysisTaskEmcalJet(name, kTRUE),
106 fConeMaxPhi(TMath::Pi()*2),
109 fCaloClustersCont(0),
112 fh2PtVsMassRCExLJDPhi(0),
114 fh2PtVsMassPerpConeLJ(0),
115 fpPtVsMassPerpConeLJ(0),
116 fh2PtVsMassPerpConeTJ(0),
117 fpPtVsMassPerpConeTJ(0),
119 fh2CentVsMassRCExLJ(0),
120 fh2CentVsMassPerpConeLJ(0),
121 fh2CentVsMassPerpConeTJ(0),
123 fh2MultVsMassRCExLJ(0),
124 fh2MultVsMassPerpConeLJ(0),
125 fh2MultVsMassPerpConeTJ(0)
127 // Standard constructor.
129 fh2PtVsMassRC = new TH2F*[fNcentBins];
130 fpPtVsMassRC = new TProfile*[fNcentBins];
131 fh2PtVsMassRCExLJDPhi = new TH3F*[fNcentBins];
132 fpPtVsMassRCExLJ = new TProfile*[fNcentBins];
133 fh2PtVsMassPerpConeLJ = new TH2F*[fNcentBins];
134 fpPtVsMassPerpConeLJ = new TProfile*[fNcentBins];
135 fh2PtVsMassPerpConeTJ = new TH2F*[fNcentBins];
136 fpPtVsMassPerpConeTJ = new TProfile*[fNcentBins];
138 for (Int_t i = 0; i < fNcentBins; i++) {
139 fh2PtVsMassRC[i] = 0;
141 fh2PtVsMassRCExLJDPhi[i] = 0;
142 fpPtVsMassRCExLJ[i] = 0;
143 fh2PtVsMassPerpConeLJ[i] = 0;
144 fpPtVsMassPerpConeLJ[i] = 0;
145 fh2PtVsMassPerpConeTJ[i] = 0;
146 fpPtVsMassPerpConeTJ[i] = 0;
149 SetMakeGeneralHistograms(kTRUE);
152 //________________________________________________________________________
153 AliAnalysisTaskEmcalJetMassBkg::~AliAnalysisTaskEmcalJetMassBkg()
158 //________________________________________________________________________
159 void AliAnalysisTaskEmcalJetMassBkg::UserCreateOutputObjects()
161 // Create user output.
163 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
165 fJetsCont = GetJetContainer(fContainerBase);
166 fTracksCont = fJetsCont->GetParticleContainer();
167 fCaloClustersCont = fJetsCont->GetClusterContainer();
169 Bool_t oldStatus = TH1::AddDirectoryStatus();
170 TH1::AddDirectory(kFALSE);
172 const Int_t nBinsPt = 250;
173 const Double_t minPt = -50.;
174 const Double_t maxPt = 200.;
176 const Int_t nBinsCent = 100;
177 const Double_t minCent = 0.;
178 const Double_t maxCent = 100.;
180 const Int_t nBinsMult = 400;
181 const Double_t minMult = 0.;
182 const Double_t maxMult = 4000.;
184 fh2CentVsMassRC = new TH2F("fh2CentVsMassRC","fh2CentVsMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
185 fOutput->Add(fh2CentVsMassRC);
187 fh2CentVsMassRCExLJ = new TH2F("fh2CentVsMassRCExLJ","fh2CentVsMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
188 fOutput->Add(fh2CentVsMassRCExLJ);
190 fh2CentVsMassPerpConeLJ = new TH2F("fh2CentVsMassPerpConeLJ","fh2CentVsMassPerpConeLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
191 fOutput->Add(fh2CentVsMassPerpConeLJ);
193 fh2CentVsMassPerpConeTJ = new TH2F("fh2CentVsMassPerpConeTJ","fh2CentVsMassPerpConeTJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
194 fOutput->Add(fh2CentVsMassPerpConeTJ);
196 fh2MultVsMassRC = new TH2F("fh2MultVsMassRC","fh2MultVsMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
197 fOutput->Add(fh2MultVsMassRC);
199 fh2MultVsMassRCExLJ = new TH2F("fh2MultVsMassRCExLJ","fh2MultVsMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
200 fOutput->Add(fh2MultVsMassRCExLJ);
202 fh2MultVsMassPerpConeLJ = new TH2F("fh2MultVsMassPerpConeLJ","fh2MultVsMassPerpConeLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
203 fOutput->Add(fh2MultVsMassPerpConeLJ);
205 fh2MultVsMassPerpConeTJ = new TH2F("fh2MultVsMassPerpConeTJ","fh2MultVsMassPerpConeTJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
206 fOutput->Add(fh2MultVsMassPerpConeTJ);
208 TString histName = "";
209 TString histTitle = "";
210 for (Int_t i = 0; i < fNcentBins; i++) {
211 histName = TString::Format("fh2PtVsMassRC_%d",i);
212 histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
213 fh2PtVsMassRC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
214 fOutput->Add(fh2PtVsMassRC[i]);
216 histName = TString::Format("fpPtVsMassRC_%d",i);
217 histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
218 fpPtVsMassRC[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
219 fOutput->Add(fpPtVsMassRC[i]);
221 histName = TString::Format("fh2PtVsMassRCExLJDPhi_%d",i);
222 histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
223 fh2PtVsMassRCExLJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,72,-0.5*TMath::Pi(),1.5*TMath::Pi());
224 fOutput->Add(fh2PtVsMassRCExLJDPhi[i]);
226 histName = TString::Format("fpPtVsMassRCExLJ_%d",i);
227 histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
228 fpPtVsMassRCExLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
229 fOutput->Add(fpPtVsMassRCExLJ[i]);
231 histName = TString::Format("fh2PtVsMassPerpConeLJ_%d",i);
232 histTitle = TString::Format("%s;#it{p}_{T,PerpConeLJ};#it{M}_{PerpConeLJ}",histName.Data());
233 fh2PtVsMassPerpConeLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
234 fOutput->Add(fh2PtVsMassPerpConeLJ[i]);
236 histName = TString::Format("fpPtVsMassPerpConeLJ_%d",i);
237 histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
238 fpPtVsMassPerpConeLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
239 fOutput->Add(fpPtVsMassPerpConeLJ[i]);
241 histName = TString::Format("fh2PtVsMassPerpConeTJ_%d",i);
242 histTitle = TString::Format("%s;#it{p}_{T,PerpConeTJ};#it{M}_{PerpConeTJ}",histName.Data());
243 fh2PtVsMassPerpConeTJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
244 fOutput->Add(fh2PtVsMassPerpConeTJ[i]);
246 histName = TString::Format("fpPtVsMassPerpConeTJ_%d",i);
247 histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
248 fpPtVsMassPerpConeTJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
249 fOutput->Add(fpPtVsMassPerpConeTJ[i]);
252 // =========== Switch on Sumw2 for all histos ===========
253 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
254 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
259 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
263 TH1::AddDirectory(oldStatus);
265 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
268 //________________________________________________________________________
269 Bool_t AliAnalysisTaskEmcalJetMassBkg::Run()
271 // Run analysis code here, if needed. It will be executed before FillHistograms().
276 //________________________________________________________________________
277 Bool_t AliAnalysisTaskEmcalJetMassBkg::FillHistograms()
281 const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
282 Double_t rho = GetRhoVal(fContainerBase);
283 Int_t trackMult = fTracksCont->GetNAcceptedParticles();
285 TLorentzVector lvRC(0.,0.,0.,0.);
290 for (Int_t i = 0; i < fRCperEvent; i++) {
291 // Simple random cones
292 lvRC.SetPxPyPzE(0.,0.,0.,0.);
296 GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
299 fh2PtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
300 fpPtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
301 fh2CentVsMassRC->Fill(fCent,RCmass);
302 fh2MultVsMassRC->Fill(trackMult,RCmass);
307 // Random cones far from leading jet(s)
308 AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
309 lvRC.SetPxPyPzE(0.,0.,0.,0.);
313 GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
315 if (RCpt > 0 && jet) {
316 Float_t dphi = RCphi - jet->Phi();
317 if (dphi > 1.5*TMath::Pi()) dphi -= TMath::Pi() * 2;
318 if (dphi < -0.5*TMath::Pi()) dphi += TMath::Pi() * 2;
319 fh2PtVsMassRCExLJDPhi[fCentBin]->Fill(RCpt - rho*rcArea,RCmass,dphi);
320 fpPtVsMassRCExLJ[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
321 fh2CentVsMassRCExLJ->Fill(fCent,RCmass);
322 fh2MultVsMassRCExLJ->Fill(trackMult,RCmass);
328 //cone perpendicular to leading jet
329 TLorentzVector lvPC(0.,0.,0.,0.);
334 AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
336 GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
339 fh2PtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
340 fpPtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
341 fh2CentVsMassPerpConeLJ->Fill(fCent,PCmass);
342 fh2MultVsMassPerpConeLJ->Fill(trackMult,PCmass);
345 //cone perpendicular to all tagged jets
346 for(int i = 0; i < fJetsCont->GetNJets();++i) {
347 jet = static_cast<AliEmcalJet*>(fJetsCont->GetAcceptJet(i));
350 if(jet->GetTagStatus()<1)
353 lvPC.SetPxPyPzE(0.,0.,0.,0.);
357 GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
360 fh2PtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
361 fpPtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
362 fh2CentVsMassPerpConeTJ->Fill(fCent,PCmass);
363 fh2MultVsMassPerpConeTJ->Fill(trackMult,PCmass);
372 //________________________________________________________________________
373 void AliAnalysisTaskEmcalJetMassBkg::GetRandomCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi,
374 AliParticleContainer* tracks, AliClusterContainer* clusters,
375 AliEmcalJet *jet) const
378 lvRC.SetPxPyPzE(0.,0.,0.,0.);
384 if (!tracks && !clusters)
395 Float_t maxEta = fConeMaxEta;
396 Float_t minEta = fConeMinEta;
397 Float_t maxPhi = fConeMaxPhi;
398 Float_t minPhi = fConeMinPhi;
400 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
401 if (minPhi < 0) minPhi = 0;
405 Bool_t reject = kTRUE;
407 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
408 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
409 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
412 } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
414 if (repeats == 999) {
415 AliWarning(Form("%s: Could not get random cone!", GetName()));
419 GetCone(lvRC,pt,eta,phi,tracks,clusters);
424 //________________________________________________________________________
425 void AliAnalysisTaskEmcalJetMassBkg::GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const
429 lvRC.SetPxPyPzE(0.,0.,0.,0.);
432 AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
434 TLorentzVector nPart;
435 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
437 Float_t cluseta = nPart.Eta();
438 Float_t clusphi = nPart.Phi();
440 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
441 clusphi += 2 * TMath::Pi();
442 if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
443 clusphi -= 2 * TMath::Pi();
445 Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
446 if (d <= fConeRadius) {
448 TLorentzVector lvcl(nPart.Px(),nPart.Py(),nPart.Pz(),nPart.E());
452 cluster = clusters->GetNextAcceptCluster();
457 AliVParticle* track = tracks->GetNextAcceptParticle(0);
459 Float_t tracketa = track->Eta();
460 Float_t trackphi = track->Phi();
462 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
463 trackphi += 2 * TMath::Pi();
464 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
465 trackphi -= 2 * TMath::Pi();
467 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
468 if (d <= fConeRadius) {
470 TLorentzVector lvtr(track->Px(),track->Py(),track->Pz(),track->E());
474 track = tracks->GetNextAcceptParticle();
480 //________________________________________________________________________
481 void AliAnalysisTaskEmcalJetMassBkg::GetPerpCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet) const
484 lvRC.SetPxPyPzE(0.,0.,0.,0.);
490 if (!tracks && !clusters)
496 Float_t LJeta = jet->Eta();
497 Float_t LJphi = jet->Phi();
500 phi = LJphi + 0.5*TMath::Pi();
501 if(phi>TMath::TwoPi()) phi-=TMath::TwoPi();
503 GetCone(lvRC,pt,eta,phi,tracks,clusters);
506 //________________________________________________________________________
507 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiEMCAL()
509 // Set default cuts for full cones
511 SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
512 SetConePhiLimits(1.405+fConeRadius,3.135-fConeRadius);
515 //________________________________________________________________________
516 void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiTPC()
518 // Set default cuts for charged cones
520 SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
521 SetConePhiLimits(-10, 10);
524 //________________________________________________________________________
525 void AliAnalysisTaskEmcalJetMassBkg::ExecOnce() {
527 AliAnalysisTaskEmcalJet::ExecOnce();
529 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
530 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
532 if (fRCperEvent < 0) {
533 Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
534 Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
535 fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
536 if (fRCperEvent == 0)
541 fMinRC2LJ = fConeRadius * 1.5;
543 const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
544 if (fMinRC2LJ > maxDist) {
545 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
546 "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
552 //________________________________________________________________________
553 Bool_t AliAnalysisTaskEmcalJetMassBkg::RetrieveEventObjects() {
555 // retrieve event objects
558 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
565 //_______________________________________________________________________
566 void AliAnalysisTaskEmcalJetMassBkg::Terminate(Option_t *)
568 // Called once at the end of the analysis.