AliPhysics  35e5fca (35e5fca)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetHF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // general ROOT includes
17 #include <TCanvas.h>
18 #include <TChain.h>
19 #include <TMath.h>
20 #include <TProfile.h>
21 #include <TAxis.h>
22 #include <TClonesArray.h>
23 #include <TH1F.h>
24 #include <TH2F.h>
25 #include <TH3F.h>
26 #include <THnSparse.h>
27 #include <TList.h>
28 #include <TLorentzVector.h>
29 #include <TParameter.h>
30 #include <TParticle.h>
31 #include <TTree.h>
32 #include <TVector3.h>
33 #include <TObjArray.h>
34 
35 #include <AliVCluster.h>
36 #include <AliVParticle.h>
37 #include <AliLog.h>
38 
39 #include "AliAnalysisTask.h"
40 #include "AliAnalysisManager.h"
41 
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliESDtrackCuts.h"
45 #include "AliAODEvent.h"
46 #include "AliAODHandler.h"
47 
48 #include "AliTLorentzVector.h"
49 #include "AliEmcalJet.h"
50 #include "AliRhoParameter.h"
51 #include "AliJetContainer.h"
52 #include "AliParticleContainer.h"
53 #include "AliClusterContainer.h"
54 
55 #include "AliPID.h"
56 #include "AliESDpid.h"
57 #include "AliAODPid.h"
58 #include "AliPIDResponse.h"
59 #include "AliKFParticle.h"
60 #include "AliKFVertex.h"
61 
63 
64 using std::cout;
65 using std::endl;
66 
67 
71 
77 fVevent(),
78 fESD(),
79 fAOD(),
80 fpidResponse(),
81 fEventCounter(),
82 fInvmassCut(),
83 fHistManager(),
84 fdEdx(), //Histo's
85 fM20(),
86 fM02(),
87 fM20EovP(),
88 fM02EovP(),
89 fInvmassLS(),
90 fInvmassULS(),
91 fEMCTrketa(),
92 fEMCTrkphi(),
93 fHistJetEovP(),
94 fHistJetEovPvPt(),
95 fHistJetEovPvEP(),
96 fHistClusEovP(),
97 fHistClusEovPnonlin(),
98 fHistClusEovPHadCorr(),
99 fHistClusTrackMatchdPhi(),
100 fHistClusTrackMatchdEta()
101 {
102 }
103 
110 AliAnalysisTaskEmcalJet(name, kTRUE),
111 fVevent(0),
112 fESD(0),
113 fAOD(0),
114 fpidResponse(0),
115 fEventCounter(0),
116 fInvmassCut(0),
117 fHistManager(name),
118 fdEdx(0), //Histo's
119 fM20(0),
120 fM02(0),
121 fM20EovP(0),
122 fM02EovP(0),
123 fInvmassLS(0),
124 fInvmassULS(0),
125 fEMCTrketa(0),
126 fEMCTrkphi(0),
127 fHistJetEovP(0),
128 fHistJetEovPvPt(0),
129 fHistJetEovPvEP(0),
130 fHistClusEovP(0),
131 fHistClusEovPnonlin(0),
132 fHistClusEovPHadCorr(0),
133 fHistClusTrackMatchdPhi(0),
134 fHistClusTrackMatchdEta(0)
135 {
137 }
138 
143 {
144 }
145 
151 {
153 
155  //Automatic determination of the analysis mode//
157  AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
158  if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
159  SetAODAnalysis();
160  } else {
161  SetESDAnalysis();
162  }
163  printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
164 
165 
170 
171  TIter next(fHistManager.GetListOfHistograms());
172  TObject* obj = 0;
173  while ((obj = next())) {
174  fOutput->Add(obj);
175  }
176  fOutput->Add(fdEdx);
177  fOutput->Add(fM20);
178  fOutput->Add(fM02);
179  fOutput->Add(fM20EovP);
180  fOutput->Add(fM02EovP);
181  fOutput->Add(fInvmassLS);
182  fOutput->Add(fInvmassULS);
183  fOutput->Add(fEMCTrketa);
184  fOutput->Add(fEMCTrkphi);
185  fOutput->Add(fHistJetEovP);
186  fOutput->Add(fHistJetEovPvPt);
187  fOutput->Add(fHistJetEovPvEP);
188  fOutput->Add(fHistClusEovP);
193 
194 
195  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
196 }
197 
198 /*
199  * This function allocates the histograms for basic EMCal cluster QA.
200  * A set of histograms (energy, eta, phi, number of cluster) is allocated
201  * per each cluster container and per each centrality bin.
202  */
204 {
205  TString histname;
206  TString histtitle;
207  TString groupname;
208  AliClusterContainer* clusCont = 0;
209  TIter next(&fClusterCollArray);
210  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
211  groupname = clusCont->GetName();
212  fHistManager.CreateHistoGroup(groupname);
213 
214 
215  fInvmassULS = new TH1F("InvmassULS", "Inv mass of ULS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
216  fInvmassLS = new TH1F("InvmassLS", "Inv mass of LS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
217  fHistClusEovP = new TH1F("ClusEovP","Cluster EovP",200,0,1.3);
218  fHistClusEovPnonlin = new TH1F("ClusEovPnonlin","nonlin Cluster EovP",200,0,5);
219  fHistClusEovPHadCorr = new TH1F("ClusEovPHaddCorr","HadCorr Cluster EovP",200,0,5);
220  fHistClusTrackMatchdPhi = new TH1F("HistClusTrackMatchdPhi","Cluster_{Phi} - Track_{Phi}; dPhi; counts;",1000,-2.0,2.0);
221  fHistClusTrackMatchdEta = new TH1F("HistClusTrackMatchdEta","Cluster_{Eta} - Track_{Eta}; dPhi; counts;",1000,-2.0,2.0);
222 
223  for (Int_t cent = 0; cent < fNcentBins; cent++) {
224  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
225  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
226  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
227 
228  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
229  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
230  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
231 
232  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
233  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
234  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
235 
236  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
237  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
238  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
239 
240  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
241  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
242  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
243 
244  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
245  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
246  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
247 
248  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
249  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
250  if (fForceBeamType != kpp) {
251  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
252  }
253  else {
254  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
255  }
256  }
257  }
258 }
259 
260 /*
261  * This function allocates the histograms for basic EMCal QA.
262  * One 2D histogram with the cell energy spectra and the number of cells
263  * per event is allocated per each centrality bin.
264  */
266 {
267  TString histname;
268  TString histtitle;
269  TString groupname(fCaloCellsName);
270 
271  fHistManager.CreateHistoGroup(groupname);
272  for (Int_t cent = 0; cent < fNcentBins; cent++) {
273  histname = TString::Format("%s/histCellEnergyvsAbsId_%d", groupname.Data(), cent);
274  histtitle = TString::Format("%s;cell abs. ID;#it{E}_{cell} (GeV);counts", histname.Data());
275  fHistManager.CreateTH2(histname, histtitle, 20000, 0, 20000, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
276 
277  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
278  histtitle = TString::Format("%s;number of cells;events", histname.Data());
279  if (fForceBeamType != kpp) {
280  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
281  }
282  else {
283  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
284  }
285  }
286 }
287 
288 /*
289  * This function allocates the histograms for basic tracking QA.
290  * A set of histograms (pT, eta, phi, difference between kinematic properties
291  * at the vertex and at the EMCal surface, number of tracks) is allocated
292  * per each particle container and per each centrality bin.
293  */
295 {
296  TString histname;
297  TString histtitle;
298  TString groupname;
299  AliParticleContainer* partCont = 0;
300  TIter next(&fParticleCollArray);
301  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
302  groupname = partCont->GetName();
303  fHistManager.CreateHistoGroup(groupname);
304  for (Int_t cent = 0; cent < fNcentBins; cent++) {
305  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
306  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
307  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
308 
309  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
310  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
311  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
312 
313  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
314  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
315  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
316 
317  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
318  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
319  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
320  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
321 
322  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
323  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
324  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
325 
326  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
327  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{p}_{T,track}^{vertex} - #it{p}_{T,track}^{EMCal} (GeV/#it{c});counts", histname.Data());
328  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
329 
330  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
331  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
332  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
333  }
334 
335  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
336  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
337  if (fForceBeamType != kpp) {
338  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
339  }
340  else {
341  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
342  }
343  }
344  }
345 }
346 
347 /*
348  * This function allocates the histograms for basic jet QA.
349  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
350  * per each jet container and per each centrality bin.
351  */
353 {
354  TString histname;
355  TString histtitle;
356  TString groupname;
357  AliJetContainer* jetCont = 0;
358  TIter next(&fJetCollArray);
359  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
360  groupname = jetCont->GetName();
361  fHistManager.CreateHistoGroup(groupname);
362 
363  fHistJetEovP = new TH1F("JetEovP","HFE Jet EovP",200,0,5);
364  fHistJetEovPvPt = new TH2F("JetEovPvPt","HFE Jet EovP vs Pt",100,0,1,100,0,100);
365 
366  for (Int_t cent = 0; cent < fNcentBins; cent++) {
367  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
368  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
369  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
370 
371  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
372  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
373  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
374 
375  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
376  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
377  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
378 
379  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
380  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
381  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
382 
383  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
384  histtitle = TString::Format("%s;number of jets;events", histname.Data());
385  if (fForceBeamType != kpp) {
386  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
387  }
388  else {
389  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
390  }
391 
392  if (!jetCont->GetRhoName().IsNull()) {
393  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
394  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
395  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
396  }
397  }
398  }
399 }
400 
408 {
409  DoJetLoop();
410  DoTrackLoop();
411  DoClusterLoop();
412  DoCellLoop();
413 
414  return kTRUE;
415 }
416 
422 {
423  TString histname;
424  TString groupname;
425  AliJetContainer* jetCont = 0;
426  TIter next(&fJetCollArray);
427 
428  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
429  groupname = jetCont->GetName();
430  UInt_t count = 0;
431  Double_t rhoVal = 0;
432  rhoVal = jetCont->GetRhoVal();
433 
434  for(auto jet : jetCont->accepted()) {
435  if (!jet) continue;
436  count++;
437 
438  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
439  fHistManager.FillTH1(histname, jet->Pt());
440 
441  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
442  fHistManager.FillTH1(histname, jet->Area());
443 
444  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
445  fHistManager.FillTH1(histname, jet->Phi());
446 
447  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
448  fHistManager.FillTH1(histname, jet->Eta());
449 
450  if (jetCont->GetRhoParameter()) {
451  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
452  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
453  }
454  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
455  fHistManager.FillTH1(histname, count);
456 
457 
458 
459  Float_t ptLeading = jetCont->GetLeadingHadronPt(jet);
460  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
461  TLorentzVector leadPart;
462  jetCont->GetLeadingHadronMomentum(leadPart, jet);
463 
464  Double_t JetPhi = jet->Phi();
465  Double_t JetEta = jet->Eta();
466  Double_t ep = jet->Phi() - fEPV0;
467  while (ep < 0) ep += TMath::Pi();
468  while (ep >= TMath::Pi()) ep -= TMath::Pi();//Get Event Plane info
469 
470 
471  //cout<<"JET Corr Pt: "<<corrPt<<" Pt of Leading Hadron: "<< ptLeading<<endl;
472 
473  AliParticleContainer* tracks = jetCont->GetParticleContainer();
474  if (tracks) {
475  for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
476  AliVParticle *track = jet->TrackAt(it, tracks->GetArray());
477 
478  if (track) {
479  AliVTrack *Vtrack = dynamic_cast<AliVTrack*>(track);
480  //Track Quality Cuts
481  /*
482  if(fAOD)
483  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue; //mimimum cuts
484 
485  //reject kink
486  Bool_t kinkmotherpass = kTRUE;
487  for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
488  if(Vtrack->GetID() == listofmotherkink[kinkmother]) {
489  kinkmotherpass = kFALSE;
490  continue;
491  }
492  }
493  if(!kinkmotherpass) continue;
494 
495  //other cuts
496  Double_t d0z0[2]={-999,-999}, cov[3];
497  Double_t DCAxyCut = 2.4, DCAzCut = 3.2;
498  if(fAOD){
499  if(Vtrack->GetNumberOfTPCClusters() < 80) continue;
500  if(Vtrack->GetNumberOfITSClusters() < 3) continue;
501  if((!(Vtrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(Vtrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
502  if(!(Vtrack->HasPointOnITSLayer(0) || Vtrack->HasPointOnITSLayer(1))) continue;
503 
504  if(Vtrack->PropagateToDCA(pVtx, fVevent->GetMagneticField(), 20., d0z0, cov))
505  if(TMath::Abs(d0z0[0]) > DCAxyCut || TMath::Abs(d0z0[1]) > DCAzCut) continue;
506  //To be done : Add cuts to apply Chi2PerITSCls < 6 and N shared Cls ITS < 4
507  }
508  */
509 
510  Int_t JetTrackLabel = Vtrack->GetID();
511  Double_t JetTrackP = track->P();
512  Double_t JetTrackPt = track->Pt();
513  Double_t dphi = TVector2::Phi_0_2pi(track->Phi() - jet->Phi());
514  Double_t deta = track->Eta() - jet->Eta();
515  Double_t dist = TMath::Sqrt(deta * deta + dphi * dphi);
516 
517  //+++++----Track matching to EMCAL
518  Int_t EMCalIndex = -1;
519  EMCalIndex = Vtrack->GetEMCALcluster();
520  if(EMCalIndex < 0) continue;
521 
522  //AliVCluster *clustMatch = (AliVCluster*)fVevent->GetCaloCluster(EMCalIndex);
523  AliClusterContainer* ClusterMatchCont = jetCont->GetClusterContainer();
524 
525  if (!ClusterMatchCont) continue;
526  AliVCluster* clusmatch = ClusterMatchCont->GetCluster(EMCalIndex);
527 
528  if(clusmatch->E() < 0.500) continue;//Bad Cluster matching below 500 MeV
529  if(!(clusmatch && clusmatch->IsEMCAL())) continue;
530  Double_t clustMatchE = clusmatch->E();
531  Double_t clustMatchNonLinE = clusmatch->GetNonLinCorrEnergy();
532  Double_t clustMatchHadCorr = clusmatch->GetHadCorrEnergy();
533  //++++---EMCal Cluster Position----
534  TVector3 clustpos;
535  Float_t emcx[3]; // cluster pos
536  Double_t detaTrckClus = -999. , dphiTrckClus = -999. , emcphi = -999., emceta = -999. ;
537  clusmatch->GetPosition(emcx);
538  clustpos.SetXYZ(emcx[0],emcx[1],emcx[2]);
539  emcphi = clustpos.Phi();
540  emceta = clustpos.Eta();
541  Double_t TrackEtaEMC = -999., TrackPhiEMC = -999., TrackPEMC = -999., TrackPtEMC = -999.;
542 
543  //++++---Jet Electrons Finder----
544  Double_t fTPCnSigma = 0.0;
545  fTPCnSigma = fpidResponse->NumberOfSigmasTPC(Vtrack, AliPID::kElectron);
546 
547  Double_t JetElectronEovP = clustMatchNonLinE / track->P();
548  TrackEtaEMC = Vtrack->GetTrackEtaOnEMCal();
549  TrackPhiEMC = Vtrack->GetTrackPhiOnEMCal();
550  TrackPEMC = Vtrack->GetTrackPOnEMCal();
551  TrackPtEMC = Vtrack->GetTrackPtOnEMCal();
552  detaTrckClus = TrackEtaEMC - track->Eta();
553  dphiTrckClus = TrackPhiEMC - track->Phi();
554  Double_t M20 = clusmatch->GetM20();
555  Double_t M02 = clusmatch->GetM02();
556 
557  //if(fTPCnSigma > -1 && fTPCnSigma < 3 && JetElectronEovP>0.8 && JetElectronEovP<1.2 && M02 > 0.006 && M02 < 0.35)
558  //{
559  //++++---Photonic Electrons----
560  Bool_t fFlagPhotonicElec = kFALSE;
561 
562  Double_t Me = 0.0005109989;
563  Double_t Me2 = Me * Me;
564 
565  AliParticleContainer* partContHFE = 0;
566  TIter next(&fParticleCollArray);
567 
568  SelectPhotonicElectron(JetTrackLabel, Vtrack, fFlagPhotonicElec);
569 
570  /*
571  while ((partContHFE = static_cast<AliParticleContainer*>(next()))) {
572  for(auto part : partContHFE->accepted()) {
573  if (!part) continue;
574  if (partContHFE->GetLoadedClass()->InheritsFrom("AliVTrack")) {
575  const AliVTrack* trackPhotonic = static_cast<const AliVTrack*>(part);
576 
577  SelectPhotonicElectron(trackPhotonic,Vtrack,fFlagPhotonicElec);
578  }
579  }//Loop over all particles in Event to look for photonic track match to Jet electron
580  }
581  */
582 
583  //++++---Jet HFE Electron Correlations----
584  Double_t dPhiJetTrk = -999., dEtaJetPhi = -999.;
585 
586 
587  fHistJetEovP->Fill(clustMatchNonLinE / Vtrack->P());
588 
589 
590  fHistClusTrackMatchdEta->Fill(detaTrckClus);
591  fHistClusTrackMatchdPhi->Fill(emcphi - TrackPhiEMC);
592 
593  //++++---electron selection cuts----
594  //if(fTPCnSigma > -1 && fTPCnSigma < 3 && eop>0.8 && eop<1.2 && m02 > 0.006 && m02 < 0.35)
595  //}
596 
597  }// Accepted Jet Track
598  }//Track loop
599  }//Track Container
600 
601 
602  //Jet Clusters
603  /*
604  AliClusterContainer* clusters = jets->GetClusterContainer();
605  if (clusters) {
606  for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
607  AliVCluster *cluster = jet->ClusterAt(ic, clusters->GetArray());
608 
609  if (cluster) {
610  TLorentzVector nPart;
611  if (clusters->GetDefaultClusterEnergy() >=0 && clusters->GetDefaultClusterEnergy() < AliVCluster::kLastUserDefEnergy) {
612  cluster->GetMomentum(nPart, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
613  }
614  else {
615  cluster->GetMomentum(nPart, fVertex);
616  }
617 
618 
619  }
620  }//loop over clusters
621  }//Get Cluster Container
622  */
623 
624 
625  }//jet loop
626  }//Jet Container
627 }//DoJetLoop()
628 
634 {
636 
637  TString histname;
638  TString groupname;
639  AliParticleContainer* partCont = 0;
640  TIter next(&fParticleCollArray);
641  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
642  groupname = partCont->GetName();
643  UInt_t count = 0;
644  for(auto part : partCont->accepted()) {
645  if (!part) continue;
646  count++;
647 
648  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
649  fHistManager.FillTH1(histname, part->Pt());
650 
651  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
652  fHistManager.FillTH1(histname, part->Phi());
653 
654  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
655  fHistManager.FillTH1(histname, part->Eta());
656 
657  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
658  const AliVTrack* track = static_cast<const AliVTrack*>(part);
659 
660  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
661  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
662 
663  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
664  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
665 
666  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
667  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
668 
669  Double_t TrackP = track->P();
670 
671  //----+++Identify Non-HFE
672  //SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
673 
674  if (clusCont) {
675  Int_t iCluster = track->GetEMCALcluster();
676  if (iCluster >= 0) {
677  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
678  if (cluster) {
679 
680  Double_t ClusterE = cluster->E();
681  Double_t ClusterNonLinE = cluster->GetNonLinCorrEnergy();
682  Double_t ClusterHadCorr = cluster->GetHadCorrEnergy();
683 
684  Double_t sigma = 0.0;
685  sigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
686  fHistClusEovP->Fill(ClusterE / TrackP);
687  fHistClusEovPnonlin->Fill( ClusterNonLinE / TrackP);
688  fHistClusEovPHadCorr->Fill( ClusterHadCorr / TrackP);
689 
690 
691 
692  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
693  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
694  }
695  }
696  }
697  }
698  }
699 
700  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
701  fHistManager.FillTH1(histname, count);
702  }
703 }
704 
710 {
711  TString histname;
712  TString groupname;
713  AliClusterContainer* clusCont = 0;
714  TIter next(&fClusterCollArray);
715  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
716  groupname = clusCont->GetName();
717 
718  for(auto cluster : clusCont->all()) {
719  if (!cluster) continue;
720 
721  if (cluster->GetIsExotic()) {
722  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
723  fHistManager.FillTH1(histname, cluster->E());
724  }
725  }
726 
727  UInt_t count = 0;
728  for(auto cluster : clusCont->accepted()) {
729  if (!cluster) continue;
730  count++;
731 
732  AliTLorentzVector nPart;
733  cluster->GetMomentum(nPart, fVertex);
734 
735  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
736  fHistManager.FillTH1(histname, cluster->E());
737 
738  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
739  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
740 
741  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
742  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
743 
744  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
745  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
746 
747  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
748  fHistManager.FillTH1(histname, nPart.Eta());
749  }
750 
751  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
752  fHistManager.FillTH1(histname, count);
753  }
754 }
755 
761 {
762  if (!fCaloCells) return;
763 
764  TString histname;
765 
766  const Short_t ncells = fCaloCells->GetNumberOfCells();
767 
768  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
769  fHistManager.FillTH1(histname, ncells);
770 
771  histname = TString::Format("%s/histCellEnergyvsAbsId_%d", fCaloCellsName.Data(), fCentBin);
772  for (Short_t pos = 0; pos < ncells; pos++) {
773  Double_t amp = fCaloCells->GetAmplitude(pos);
774  Short_t absId = fCaloCells->GetCellNumber(pos);
775 
776  fHistManager.FillTH2(histname, absId, amp);
777  }
778 }
779 
785 {
787 }
788 
797 {
798  UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
799 
800  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
801  if (!fVevent) {
802  printf("ERROR: fVEvent not available\n");
803  return kFALSE;
804  }
805 
806  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
807  if (fESD) {
808  // printf("fESD available\n");
809  //return;
810  }
811 
812  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
813  if (fAOD) {
814  // printf("fAOD available\n");
815  //return;
816  }
817 
818  const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
819 
820  fpidResponse = fInputHandler->GetPIDResponse();
821  if (!fpidResponse) {
822  printf("ERROR: fpidResponse not available\n");
823  return kFALSE;
824  }
825 
826  Int_t ntracks;
827  ntracks = fVevent->GetNumberOfTracks();
828  //printf("There are %d tracks in this event\n",ntracks);
829 
830  Double_t Zvertex = -100, Xvertex = -100, Yvertex = -100;
831  Double_t NcontV = pVtx->GetNContributors();
832 
833  // if(NcontV<2)return kTRUE; //Events with 2 Tracks
834 
836  //EMcal cluster info//
839 
840 
841  Int_t numberofvertices = 100;
842  if(fAOD) numberofvertices = fAOD->GetNumberOfVertices();
843  Double_t listofmotherkink[numberofvertices];
844  Int_t numberofmotherkink = 0;
845  if(IsAODanalysis())
846  {
847  for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
848  AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
849  if(!aodvertex) continue;
850  if(aodvertex->GetType()==AliAODVertex::kKink) {
851  AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
852  if(!mother) continue;
853  Int_t idmother = mother->GetID();
854  listofmotherkink[numberofmotherkink] = idmother;
855  numberofmotherkink++;
856  }
857  }
858  } //+++
859 
860  fEventCounter++;
861 
862  return kTRUE;
863 }
864 
865 //________________________________________________________________________
867 {
869  //EMCAL cluster information//
871 
872  Int_t Nclust = -999;
873  TVector3 clustpos;
874  Float_t emcx[3]; // cluster pos
875  Double_t clustE=-999, emcphi = -999, emceta=-999;
876  Nclust = fVevent->GetNumberOfCaloClusters();
877  for(Int_t icl=0; icl<Nclust; icl++)
878  {
879  AliVCluster *clust = 0x0;
880  clust = fVevent->GetCaloCluster(icl);
881  if(!clust) printf("ERROR: Could not receive cluster matched calibrated from track %d\n", icl);
882 
883  if(clust && clust->IsEMCAL())
884  {
885  clustE = clust->E();
886  clust->GetPosition(emcx);
887  clustpos.SetXYZ(emcx[0],emcx[1],emcx[2]);
888  emcphi = clustpos.Phi();
889  emceta = clustpos.Eta();
890 
891  }
892  }
893 }
894 //________________________________________________________________________
895 
896 //___________________________________________________________________________
897 
898 void AliAnalysisTaskEmcalJetHF::SelectPhotonicElectron(Int_t itrack, AliVTrack *track, Bool_t &fFlagPhotonicElec)
899 {
901  //Photonic electron selection//
903 
904  Bool_t flagPhotonicElec = kFALSE;
905  Double_t ptAsso=-999., nsigma=-999.0;
906  Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
907 
908  for(Int_t jTracks = 0; jTracks<fVevent->GetNumberOfTracks(); jTracks++){
909 
910  AliVParticle* VtrackAsso = fVevent->GetTrack(jTracks);
911  if (!VtrackAsso) {
912  printf("ERROR: Could not receive track %d\n", jTracks);
913  continue;
914  }
915  AliVTrack *trackAsso = dynamic_cast<AliVTrack*>(VtrackAsso);
916  if(!trackAsso) continue;
917  if(trackAsso->GetID() == itrack){
918  //cout << "Matched track to itself: " << trackAsso->GetID() <<" " << itrack<<endl;
919  continue;
920  }
921 
922  if(IsAODanalysis()) {
923  AliAODTrack *atrackAsso = dynamic_cast<AliAODTrack*>(VtrackAsso);
924  if(!atrackAsso) continue;
925  if(!atrackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
926  if(atrackAsso->GetTPCNcls() < 70) continue;
927  if(!(atrackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
928  if(!(atrackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
929  }
930 
931  nsigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
932  ptAsso = trackAsso->Pt();
933  Int_t chargeAsso = trackAsso->Charge();
934  Int_t charge = track->Charge();
935 
936  if(ptAsso <0.2) continue;
937  if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
938  if(nsigma < -3 || nsigma > 3) continue;
939 
940  Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
941  if(charge>0) fPDGe1 = -11;
942  if(chargeAsso>0) fPDGe2 = -11;
943 
944  if(charge == chargeAsso) fFlagLS = kTRUE;
945  if(charge != chargeAsso) fFlagULS = kTRUE;
946 
947  AliKFParticle::SetField(fVevent->GetMagneticField());
948 
949  AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
950  AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
951  AliKFParticle recg(ge1, ge2);
952 
953  if(recg.GetNDF()<1) continue;
954  Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
955  if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
956 
957  Double_t mass=-999., width = -999;
958  Int_t MassCorrect;
959  MassCorrect = recg.GetMass(mass,width);
960 
961  if(fFlagLS && track->Pt()>1) fInvmassLS->Fill(mass);
962  if(fFlagULS && track->Pt()>1) fInvmassULS->Fill(mass);
963 
964  if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
965  flagPhotonicElec = kTRUE;
966  }
967  }
968  fFlagPhotonicElec = flagPhotonicElec;
969 }
970 
971 
976 {
977  cout<<"#######################"<<endl;
978  cout<<"#### Task Finished ####"<<endl;
979  cout<<"#######################"<<endl;
980 }
Int_t charge
TH1F * fHistJetEovPvEP
Jet EovP vs Pt.
TH1F * fHistClusTrackMatchdPhi
EovP using Hadron Corected Energy.
THashList * CreateHistoGroup(const char *groupname)
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
AliAODEvent * fAOD
ESD object.
Double_t mass
TH2F * fdEdx
Invariant mass Cut for photons.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Declaration of class AliTLorentzVector.
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
AliClusterContainer * GetClusterContainer() const
Int_t fCentBin
!event centrality bin
TH1F * fHistClusEovPHadCorr
EovP using Non linear corrected energy.
void SelectPhotonicElectron(Int_t itrack, AliVTrack *track, Bool_t &fFlagPhotonicElec)
THistManager fHistManager
Histogram manager.
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
particle/track collection array
const AliClusterIterableContainer all() const
Double_t * sigma
AliParticleContainer * GetParticleContainer() const
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
TH1F * fHistClusTrackMatchdEta
dPhi = Phi_cluster - Phi_track Matching
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
Double_t Phi_0_2pi() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliVCluster * GetAcceptCluster(Int_t i) const
const AliClusterIterableContainer accepted() const
TString fCaloCellsName
name of calo cell collection
Double_t nsigma
Declaration of class AliAnalysisTaskEmcalJetHF. Task to perform analysis on HFE JEts.
AliVCluster * GetCluster(Int_t i) const
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
TObjArray fJetCollArray
jet collection array
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Double_t fInvmassCut
Event Counter.
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1F * fHistClusEovPnonlin
EovP using Cluster Energy.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t fVertex[3]
!event vertex
TH1F * fEMCTrketa
Invmass of ULS pairs.
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
AliESDEvent * fESD
event object
const AliParticleIterableContainer accepted() const
const char Option_t
Definition: External.C:48
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Container structure for EMCAL clusters.
TH1F * fHistClusEovP
Jet EovP vs Event Plane.
TH1F * fInvmassULS
Invmass of LS pairs.
AliPIDResponse * fpidResponse
AOD object.
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins