2 // Do subtraction for jet shapes using constituents arXiv:1403.3108
6 #include <TClonesArray.h>
10 #include <THnSparse.h>
13 #include <TLorentzVector.h>
21 #include "AliVCluster.h"
22 #include "AliVTrack.h"
23 #include "AliEmcalJet.h"
24 #include "AliRhoParameter.h"
26 #include "AliEmcalParticle.h"
27 #include "AliMCEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliGenPythiaEventHeader.h"
30 #include "AliAODMCHeader.h"
31 #include "AliAnalysisManager.h"
32 #include "AliJetContainer.h"
33 #include "AliParticleContainer.h"
35 #include "AliAnalysisTaskJetShapeConst.h"
37 ClassImp(AliAnalysisTaskJetShapeConst)
39 //________________________________________________________________________
40 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() :
41 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
45 fMinFractionShared(0),
46 fSingleTrackEmb(kFALSE),
48 fJetMassVarType(kMass),
50 fJet1Vec(new TLorentzVector()),
51 fJet2Vec(new TLorentzVector()),
52 fJetSubVec(new TLorentzVector()),
62 fh3MSubPtRawDRMatch(0x0),
64 fh3MTruePtTrueDR(0x0),
65 fh3PtTrueDeltaMDR(0x0),
66 fh3PtTrueDeltaMRelDR(0x0),
69 // Default constructor.
71 fh2MSubMatch = new TH2F*[fNcentBins];
72 fh2MSubPtRawAll = new TH2F*[fNcentBins];
73 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
74 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
75 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
76 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
77 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
78 fhnMassResponse = new THnSparse*[fNcentBins];
80 for (Int_t i = 0; i < fNcentBins; i++) {
82 fh2MSubPtRawAll[i] = 0;
83 fh3MSubPtRawDRMatch[i] = 0;
84 fh3MSubPtTrueDR[i] = 0;
85 fh3MTruePtTrueDR[i] = 0;
86 fh3PtTrueDeltaMDR[i] = 0;
87 fh3PtTrueDeltaMRelDR[i] = 0;
88 fhnMassResponse[i] = 0;
91 SetMakeGeneralHistograms(kTRUE);
92 if(fCreateTree) DefineOutput(2, TTree::Class());
95 //________________________________________________________________________
96 AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
97 AliAnalysisTaskEmcalJet(name, kTRUE),
101 fMinFractionShared(0),
102 fSingleTrackEmb(kFALSE),
104 fJetMassVarType(kMass),
106 fJet1Vec(new TLorentzVector()),
107 fJet2Vec(new TLorentzVector()),
108 fJetSubVec(new TLorentzVector()),
117 fh2MSubPtRawAll(0x0),
118 fh3MSubPtRawDRMatch(0x0),
119 fh3MSubPtTrueDR(0x0),
120 fh3MTruePtTrueDR(0x0),
121 fh3PtTrueDeltaMDR(0x0),
122 fh3PtTrueDeltaMRelDR(0x0),
125 // Standard constructor.
127 fh2MSubMatch = new TH2F*[fNcentBins];
128 fh2MSubPtRawAll = new TH2F*[fNcentBins];
129 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
130 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
131 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
132 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
133 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
134 fhnMassResponse = new THnSparse*[fNcentBins];
136 for (Int_t i = 0; i < fNcentBins; i++) {
138 fh2MSubPtRawAll[i] = 0;
139 fh3MSubPtRawDRMatch[i] = 0;
140 fh3MSubPtTrueDR[i] = 0;
141 fh3MTruePtTrueDR[i] = 0;
142 fh3PtTrueDeltaMDR[i] = 0;
143 fh3PtTrueDeltaMRelDR[i] = 0;
144 fhnMassResponse[i] = 0;
147 SetMakeGeneralHistograms(kTRUE);
148 if(fCreateTree) DefineOutput(2, TTree::Class());
151 //________________________________________________________________________
152 AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
157 //________________________________________________________________________
158 void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
160 // Create user output.
162 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
164 Bool_t oldStatus = TH1::AddDirectoryStatus();
165 TH1::AddDirectory(kFALSE);
167 const Int_t nBinsPt = 200;
168 const Double_t minPt = -50.;
169 const Double_t maxPt = 150.;
172 Double_t minM = -20.;
174 if(fJetMassVarType==kRatMPt) {
181 Double_t minDM = -50.;
182 Double_t maxDM = 50.;
183 if(fJetMassVarType==kRatMPt) {
189 const Int_t nBinsDRToLJ = 20; //distance to leading jet in Pb-Pb only event
190 const Double_t minDRToLJ = 0.;
191 const Double_t maxDRToLJ = 1.;
193 //Binning for THnSparse
194 const Int_t nBinsSparse0 = 5;
195 const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
196 const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minDRToLJ};
197 const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
199 TString histName = "";
200 TString histTitle = "";
201 TString varName = "#it{M}_{jet}";
202 if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
204 for (Int_t i = 0; i < fNcentBins; i++) {
205 histName = Form("fh2MSubMatch_%d",i);
206 histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
207 fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
208 fOutput->Add(fh2MSubMatch[i]);
210 histName = Form("fh2MSubPtRawAll_%d",i);
211 histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
212 fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
213 fOutput->Add(fh2MSubPtRawAll[i]);
215 histName = Form("fh3MSubPtRawDRMatch_%d",i);
216 histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
217 fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
218 fOutput->Add(fh3MSubPtRawDRMatch[i]);
220 histName = Form("fh3MSubPtTrueDR_%d",i);
221 histTitle = Form("fh3MSubPtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
222 fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
223 fOutput->Add(fh3MSubPtTrueDR[i]);
225 histName = Form("fh3MTruePtTrueDR_%d",i);
226 histTitle = Form("fh3MTruePtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
227 fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
228 fOutput->Add(fh3MTruePtTrueDR[i]);
230 histName = Form("fh3PtTrueDeltaMDR_%d",i);
231 histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
232 fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
233 fOutput->Add(fh3PtTrueDeltaMDR[i]);
235 histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
236 histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
237 fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
238 fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
240 histName = Form("fhnMassResponse_%d",i);
241 histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#Delta R_{LJ}",i,varName.Data(),varName.Data());
242 fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
243 fOutput->Add(fhnMassResponse[i]);
246 // =========== Switch on Sumw2 for all histos ===========
247 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
248 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
253 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
257 TH1::AddDirectory(oldStatus);
261 fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
262 fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
263 fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
264 fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
265 fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
266 fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
267 fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
268 fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
269 fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
270 fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
273 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
274 if(fCreateTree) PostData(2, fTreeJetBkg);
277 //________________________________________________________________________
278 Bool_t AliAnalysisTaskJetShapeConst::Run()
280 // Run analysis code here, if needed. It will be executed before FillHistograms().
285 //________________________________________________________________________
286 Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
290 AliEmcalJet* jet1 = NULL; //AA jet
291 AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
292 // AliEmcalJet *jet1T = NULL; //tagged AA jet
293 // AliEmcalJet *jet2T = NULL; //tagged Pythia jet
294 AliEmcalJet *jetS = NULL; //subtracted jet
295 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
296 AliJetContainer *jetContS = GetJetContainer(fContainerSub);
298 //Get leading jet in Pb-Pb event without embedded objects
299 AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
300 AliEmcalJet *jetL = NULL;
301 if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
303 AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
305 jetCont->ResetCurrentID();
306 while((jet1 = jetCont->GetNextAcceptJet())) {
310 //Get constituent subtacted version of jet
313 for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
314 //if(ifound==1) continue;
315 jetS = jetContS->GetJet(i);
316 if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
318 if(ifound==1) ilab = i;
321 if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
322 if(ifound==0) jetS = 0x0;
323 else jetS = jetContS->GetJet(ilab);
325 Double_t mjetS = jetS->M();
326 Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
327 Double_t var = mjetS;
328 if(fJetMassVarType==kRatMPt) {
329 if(ptjet1>0. || ptjet1<0.) var = mjetS/ptjet1;
333 //Fill histograms for all AA jets
334 fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
336 Double_t fraction = 0.;
338 fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
339 if(fSingleTrackEmb) {
340 AliVParticle *vp = GetEmbeddedConstituent(jet1);
342 fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
346 jet2 = jet1->ClosestJet();
347 fraction = jetCont->GetFractionSharedPt(jet1);
349 if(fMinFractionShared>0.) {
350 if(fraction>fMinFractionShared) {
351 fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
358 // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
360 //Fill histograms for matched jets
361 fh2MSubMatch[fCentBin]->Fill(var,fMatch);
363 Double_t drToLJ = -1.;
364 if(jetL) drToLJ = jet1->DeltaR(jetL);
365 fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
367 Double_t var2 = -999.;
368 if(jet2->Pt()>0. || jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
369 fh3MSubPtTrueDR[fCentBin]->Fill(var,jet2->Pt(),drToLJ);
370 fh3MTruePtTrueDR[fCentBin]->Fill(var2,jet2->Pt(),drToLJ);
371 fh3PtTrueDeltaMDR[fCentBin]->Fill(jet2->Pt(),var-var2,drToLJ);
372 if(jet2->M()>0.) fh3PtTrueDeltaMRelDR[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,drToLJ);
373 Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),drToLJ};
374 fhnMassResponse[fCentBin]->Fill(varsp);
379 fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
380 if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
381 else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
382 fArea = (Float_t)jet1->Area();
383 fAreaPhi = (Float_t)jet1->AreaPhi();
384 fAreaEta = (Float_t)jet1->AreaEta();
385 fRho = (Float_t)jetCont->GetRhoVal();
386 fRhoM = (Float_t)jetCont->GetRhoMassVal();
387 fNConst = (Int_t)jet1->GetNumberOfTracks();
397 //________________________________________________________________________
398 AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
400 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
401 AliVParticle *vp = 0x0;
402 AliVParticle *vpe = 0x0; //embedded particle
404 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
405 vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
406 if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
408 else if(vp->Pt()>vpe->Pt()) vpe = vp;
411 AliDebug(11,Form("Found %d embedded particles",nc));
416 //________________________________________________________________________
417 Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
419 // retrieve event objects
422 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
425 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
426 jetCont->LoadRhoMass(InputEvent());
431 //_______________________________________________________________________
432 void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
434 // Called once at the end of the analysis.