Charged jets (pPb): Improved trackcut analysis
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetMass.cxx
1 //
2 // Jet mass 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 "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"
31
32 #include "AliAODEvent.h"
33
34 #include "AliAnalysisTaskEmcalJetMass.h"
35
36 ClassImp(AliAnalysisTaskEmcalJetMass)
37
38 //________________________________________________________________________
39 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() : 
40   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
41   fContainerBase(0),
42   fMinFractionShared(0),
43   fJetMassAvg(0),
44   fh2PtJet1VsLeadPtAllSel(0),
45   fh2PtJet1VsLeadPtTagged(0),
46   fh2PtJet1VsLeadPtTaggedMatch(0),
47   fh2PtVsMassJet1All(0),
48   fh2PtVsMassJet1Tagged(0),
49   fh2PtVsMassJet1TaggedMatch(0),
50   fpPtVsMassJet1All(0),
51   fpPtVsMassJet1Tagged(0),
52   fpPtVsMassJet1TaggedMatch(0),
53   fh2MassVsAreaJet1All(0),
54   fh2MassVsAreaJet1Tagged(0),
55   fh2MassVsAreaJet1TaggedMatch(0),
56   fh2MassVsNConstJet1All(0),
57   fh2MassVsNConstJet1Tagged(0),
58   fh2MassVsNConstJet1TaggedMatch(0),
59   fh2EtMassOverEtRSq(0)
60 {
61   // Default constructor.
62
63   fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
64   fh2PtJet1VsLeadPtTagged        = new TH2F*[fNcentBins];
65   fh2PtJet1VsLeadPtTaggedMatch   = new TH2F*[fNcentBins];
66   fh2PtVsMassJet1All             = new TH2F*[fNcentBins];
67   fh2PtVsMassJet1Tagged          = new TH2F*[fNcentBins];
68   fh2PtVsMassJet1TaggedMatch     = new TH2F*[fNcentBins];
69   fpPtVsMassJet1All              = new TProfile*[fNcentBins];
70   fpPtVsMassJet1Tagged           = new TProfile*[fNcentBins];
71   fpPtVsMassJet1TaggedMatch      = new TProfile*[fNcentBins];
72   fh2MassVsAreaJet1All           = new TH2F*[fNcentBins];
73   fh2MassVsAreaJet1Tagged        = new TH2F*[fNcentBins];
74   fh2MassVsAreaJet1TaggedMatch   = new TH2F*[fNcentBins];
75   fh2MassVsNConstJet1All         = new TH2F*[fNcentBins];
76   fh2MassVsNConstJet1Tagged      = new TH2F*[fNcentBins];
77   fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
78   fh2EtMassOverEtRSq             = new TH2F*[fNcentBins];
79
80   for (Int_t i = 0; i < fNcentBins; i++) {
81     fh2PtJet1VsLeadPtAllSel[i]        = 0;
82     fh2PtJet1VsLeadPtTagged[i]        = 0;
83     fh2PtJet1VsLeadPtTaggedMatch[i]   = 0;
84     fh2PtVsMassJet1All[i]             = 0;
85     fh2PtVsMassJet1Tagged[i]          = 0;
86     fh2PtVsMassJet1TaggedMatch[i]     = 0;
87     fpPtVsMassJet1All[i]              = 0;
88     fpPtVsMassJet1Tagged[i]           = 0;
89     fpPtVsMassJet1TaggedMatch[i]      = 0;
90     fh2MassVsAreaJet1All[i]           = 0;
91     fh2MassVsAreaJet1Tagged[i]        = 0;
92     fh2MassVsAreaJet1TaggedMatch[i]   = 0;
93     fh2MassVsNConstJet1All[i]         = 0;
94     fh2MassVsNConstJet1Tagged[i]      = 0;
95     fh2MassVsNConstJet1TaggedMatch[i] = 0;
96     fh2EtMassOverEtRSq[i]             = 0;
97   }
98
99   SetMakeGeneralHistograms(kTRUE);
100 }
101
102 //________________________________________________________________________
103 AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) : 
104   AliAnalysisTaskEmcalJet(name, kTRUE),  
105   fContainerBase(0),
106   fMinFractionShared(0),
107   fJetMassAvg(0),
108   fh2PtJet1VsLeadPtAllSel(0),
109   fh2PtJet1VsLeadPtTagged(0),
110   fh2PtJet1VsLeadPtTaggedMatch(0),
111   fh2PtVsMassJet1All(0),
112   fh2PtVsMassJet1Tagged(0),
113   fh2PtVsMassJet1TaggedMatch(0),
114   fpPtVsMassJet1All(0),
115   fpPtVsMassJet1Tagged(0),
116   fpPtVsMassJet1TaggedMatch(0),
117   fh2MassVsAreaJet1All(0),
118   fh2MassVsAreaJet1Tagged(0),
119   fh2MassVsAreaJet1TaggedMatch(0),
120   fh2MassVsNConstJet1All(0),
121   fh2MassVsNConstJet1Tagged(0),
122   fh2MassVsNConstJet1TaggedMatch(0),
123   fh2EtMassOverEtRSq(0)
124 {
125   // Standard constructor.
126
127  fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
128   fh2PtJet1VsLeadPtTagged        = new TH2F*[fNcentBins];
129   fh2PtJet1VsLeadPtTaggedMatch   = new TH2F*[fNcentBins];
130   fh2PtVsMassJet1All             = new TH2F*[fNcentBins];
131   fh2PtVsMassJet1Tagged          = new TH2F*[fNcentBins];
132   fh2PtVsMassJet1TaggedMatch     = new TH2F*[fNcentBins];
133   fpPtVsMassJet1All              = new TProfile*[fNcentBins];
134   fpPtVsMassJet1Tagged           = new TProfile*[fNcentBins];
135   fpPtVsMassJet1TaggedMatch      = new TProfile*[fNcentBins];
136   fh2MassVsAreaJet1All           = new TH2F*[fNcentBins];
137   fh2MassVsAreaJet1Tagged        = new TH2F*[fNcentBins];
138   fh2MassVsAreaJet1TaggedMatch   = new TH2F*[fNcentBins];
139   fh2MassVsNConstJet1All         = new TH2F*[fNcentBins];
140   fh2MassVsNConstJet1Tagged      = new TH2F*[fNcentBins];
141   fh2MassVsNConstJet1TaggedMatch = new TH2F*[fNcentBins];
142   fh2EtMassOverEtRSq             = new TH2F*[fNcentBins];
143
144   for (Int_t i = 0; i < fNcentBins; i++) {
145     fh2PtJet1VsLeadPtAllSel[i]        = 0;
146     fh2PtJet1VsLeadPtTagged[i]        = 0;
147     fh2PtJet1VsLeadPtTaggedMatch[i]   = 0;
148     fh2PtVsMassJet1All[i]             = 0;
149     fh2PtVsMassJet1Tagged[i]          = 0;
150     fh2PtVsMassJet1TaggedMatch[i]     = 0;
151     fpPtVsMassJet1All[i]              = 0;
152     fpPtVsMassJet1Tagged[i]           = 0;
153     fpPtVsMassJet1TaggedMatch[i]      = 0;
154     fh2MassVsAreaJet1All[i]           = 0;
155     fh2MassVsAreaJet1Tagged[i]        = 0;
156     fh2MassVsAreaJet1TaggedMatch[i]   = 0;
157     fh2MassVsNConstJet1All[i]         = 0;
158     fh2MassVsNConstJet1Tagged[i]      = 0;
159     fh2MassVsNConstJet1TaggedMatch[i] = 0;
160     fh2EtMassOverEtRSq[i]             = 0;
161   }
162
163   SetMakeGeneralHistograms(kTRUE);
164 }
165
166 //________________________________________________________________________
167 AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
168 {
169   // Destructor.
170 }
171
172 //________________________________________________________________________
173 void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
174 {
175   // Create user output.
176
177   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
178
179   Bool_t oldStatus = TH1::AddDirectoryStatus();
180   TH1::AddDirectory(kFALSE);
181
182   const Int_t nBinsPt  = 250;
183   const Double_t minPt = -50.;
184   const Double_t maxPt = 200.;
185
186   const Int_t nBinsM  = 150;
187   const Double_t minM = -50.;
188   const Double_t maxM = 100.;
189
190   const Int_t nBinsArea = 100;
191   const Double_t minArea = 0.;
192   const Double_t maxArea = 1.;
193
194   const Int_t nBinsNConst = 100;
195   const Double_t minNConst = 0.;
196   const Double_t maxNConst = 500.;
197
198   TString histName = "";
199   TString histTitle = "";
200   for (Int_t i = 0; i < fNcentBins; i++) {
201     histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
202     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
203     fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
204     fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
205
206     histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
207     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
208     fh2PtJet1VsLeadPtTagged[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
209     fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
210
211     histName = TString::Format("fh2PtJet1VsLeadPtTaggedMatch_%d",i);
212     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
213     fh2PtJet1VsLeadPtTaggedMatch[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
214     fOutput->Add(fh2PtJet1VsLeadPtTaggedMatch[i]);
215
216     histName = TString::Format("fh2PtVsMassJet1All_%d",i);
217     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
218     fh2PtVsMassJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
219     fOutput->Add(fh2PtVsMassJet1All[i]);
220
221     histName = TString::Format("fh2PtVsMassJet1Tagged_%d",i);
222     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
223     fh2PtVsMassJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
224     fOutput->Add(fh2PtVsMassJet1Tagged[i]);
225
226     histName = TString::Format("fh2PtVsMassJet1TaggedMatch_%d",i);
227     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
228     fh2PtVsMassJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
229     fOutput->Add(fh2PtVsMassJet1TaggedMatch[i]);
230
231     histName = TString::Format("fpPtVsMassJet1All_%d",i);
232     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
233     fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
234     fOutput->Add(fpPtVsMassJet1All[i]);
235
236     histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
237     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
238     fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
239     fOutput->Add(fpPtVsMassJet1Tagged[i]);
240
241     histName = TString::Format("fpPtVsMassJet1TaggedMatch_%d",i);
242     histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
243     fpPtVsMassJet1TaggedMatch[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsM,minM,maxM);
244     fOutput->Add(fpPtVsMassJet1TaggedMatch[i]);
245
246     histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
247     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
248     fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
249     fOutput->Add(fh2MassVsAreaJet1All[i]);
250
251     histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
252     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
253     fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
254     fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
255
256     histName = TString::Format("fh2MassVsAreaJet1TaggedMatch_%d",i);
257     histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
258     fh2MassVsAreaJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsArea,minArea,maxArea);
259     fOutput->Add(fh2MassVsAreaJet1TaggedMatch[i]);
260
261     histName = TString::Format("fh2MassVsNConstJet1All_%d",i);
262     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
263     fh2MassVsNConstJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsNConst,minNConst,maxNConst);
264     fOutput->Add(fh2MassVsNConstJet1All[i]);
265
266     histName = TString::Format("fh2MassVsNConstJet1Tagged_%d",i);
267     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
268     fh2MassVsNConstJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
269     fOutput->Add(fh2MassVsNConstJet1Tagged[i]);
270
271     histName = TString::Format("fh2MassVsNConstJet1TaggedMatch_%d",i);
272     histTitle = TString::Format("%s;#it{M}_{jet1};#it{N}_{constituents}",histName.Data());
273     fh2MassVsNConstJet1TaggedMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsNConst,minNConst,maxNConst);
274     fOutput->Add(fh2MassVsNConstJet1TaggedMatch[i]);
275
276     histName = TString::Format("fh2EtMassOverEtRSq_%d",i);
277     histTitle = TString::Format("%s;#it{E}_{T};(#it{M}/(#it{E}_{T}#it{R}))^{2}",histName.Data());
278     fh2EtMassOverEtRSq[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,1.);
279     fOutput->Add(fh2EtMassOverEtRSq[i]);
280   }
281
282   // =========== Switch on Sumw2 for all histos ===========
283   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
284     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
285     if (h1){
286       h1->Sumw2();
287       continue;
288     }
289     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
290     if(hn)hn->Sumw2();
291   }
292
293   TH1::AddDirectory(oldStatus);
294
295   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
296 }
297
298 //________________________________________________________________________
299 Bool_t AliAnalysisTaskEmcalJetMass::Run()
300 {
301   // Run analysis code here, if needed. It will be executed before FillHistograms().
302
303   return kTRUE;
304 }
305
306 //________________________________________________________________________
307 Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
308 {
309   // Fill histograms.
310
311   AliEmcalJet* jet1 = NULL;
312   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
313   if(jetCont) {
314     jetCont->ResetCurrentID();
315     while((jet1 = jetCont->GetNextAcceptJet())) {
316
317       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
318       Double_t mJet1 = GetJetMass(jet1);
319       fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
320       fh2PtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
321       fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,mJet1);
322       fh2MassVsAreaJet1All[fCentBin]->Fill(mJet1,jet1->Area());
323       fh2MassVsNConstJet1All[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
324       
325       if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
326         continue;
327
328       fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
329       fh2PtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
330       fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,mJet1);
331       fh2MassVsAreaJet1Tagged[fCentBin]->Fill(mJet1,jet1->Area());
332       fh2MassVsNConstJet1Tagged[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
333
334       Double_t fraction = jetCont->GetFractionSharedPt(jet1);
335       if(fMinFractionShared>0. && fraction>fMinFractionShared) {
336         fh2PtJet1VsLeadPtTaggedMatch[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
337         fh2PtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
338         fpPtVsMassJet1TaggedMatch[fCentBin]->Fill(ptJet1,mJet1);
339         fh2MassVsAreaJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->Area());
340         fh2MassVsNConstJet1TaggedMatch[fCentBin]->Fill(mJet1,jet1->GetNumberOfConstituents());
341       }
342       
343       Double_t Et2 = jet1->M()*jet1->M() + jet1->Pt()*jet1->Pt();
344       Double_t Et = 0.;    Double_t massOverEtR = 0.;
345       if(Et2>0.) Et = TMath::Sqrt(Et2);
346       if((Et*jetCont->GetJetRadius())>0.) 
347         massOverEtR = jet1->M()/(Et*jetCont->GetJetRadius());
348       fh2EtMassOverEtRSq[fCentBin]->Fill(Et,massOverEtR*massOverEtR);
349     }
350   }
351   
352   return kTRUE;
353 }
354
355 //________________________________________________________________________
356 Double_t AliAnalysisTaskEmcalJetMass::GetJetMass(AliEmcalJet *jet) {
357   //calc subtracted jet mass
358   TLorentzVector vpS = GetSubtractedVector(jet);
359   return vpS.M();
360 }
361
362
363 //________________________________________________________________________
364 TLorentzVector AliAnalysisTaskEmcalJetMass::GetSubtractedVector(AliEmcalJet *jet) {
365   //get subtracted vector
366   TLorentzVector vpS;
367   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
368   TLorentzVector vpBkg = GetBkgVector(jet,jetCont);
369   vpS.SetPxPyPzE(jet->Px()-vpBkg.Px(),jet->Py()-vpBkg.Py(),jet->Pz()-vpBkg.Pz(),jet->E()-vpBkg.E());
370   return vpS;
371 }
372
373 //________________________________________________________________________
374 TLorentzVector AliAnalysisTaskEmcalJetMass::GetBkgVector(AliEmcalJet *jet, AliJetContainer *cont) {
375   //get background vector
376
377   Double_t rho  = cont->GetRhoVal();
378   Double_t rhom = cont->GetRhoMassVal();
379   TLorentzVector vpB;
380   Double_t aphi = jet->AreaPhi();
381   Double_t aeta = jet->AreaEta();
382   vpB.SetPxPyPzE(rho*TMath::Cos(aphi)*jet->Area(),rho*TMath::Sin(aphi)*jet->Area(),(rho+rhom)*TMath::SinH(aeta)*jet->Area(),(rho+rhom)*TMath::CosH(aeta)*jet->Area());
383   return vpB;
384 }
385
386 //________________________________________________________________________
387 Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
388   //
389   // retrieve event objects
390   //
391
392   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
393     return kFALSE;
394
395   return kTRUE;
396
397 }
398
399 //_______________________________________________________________________
400 void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *) 
401 {
402   // Called once at the end of the analysis.
403 }
404