AliPhysics  a6017e1 (a6017e1)
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 fHistManager(),
78 fVevent(),
79 fESD(),
80 fAOD(),
81 fpidResponse(),
82 fEventCounter(),
83 fInvmassCut(),
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 fHistManager(name),
112 fVevent(0),
113 fESD(0),
114 fAOD(0),
115 fpidResponse(0),
116 fEventCounter(0),
117 fInvmassCut(0),
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),
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)
Create a new group of histograms within a parent group.
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="")
Fill a 2D histogram within the container.
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="")
Create a new TH2 within the container.
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
Get the list of histograms.
Definition: THistManager.h:671
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="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
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="")
Fill a 1D histogram within the container.
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.
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
void UserCreateOutputObjects()
Main initialization function on the worker.
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