AliPhysics  c0d7b22 (c0d7b22)
 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 fHistClusEovP(),
96 fHistClusEovPnonlin(),
97 fHistClusEovPHadCorr()
98 {
99 }
100 
107 AliAnalysisTaskEmcalJet(name, kTRUE),
108 fVevent(0),
109 fESD(0),
110 fAOD(0),
111 fpidResponse(0),
112 fEventCounter(0),
113 fInvmassCut(0),
114 fHistManager(name),
115 fdEdx(0), //Histo's
116 fM20(0),
117 fM02(0),
118 fM20EovP(0),
119 fM02EovP(0),
120 fInvmassLS(0),
121 fInvmassULS(0),
122 fEMCTrketa(0),
123 fEMCTrkphi(0),
124 fHistJetEovP(0),
125 fHistJetEovPvPt(0),
126 fHistClusEovP(0),
127 fHistClusEovPnonlin(0),
128 fHistClusEovPHadCorr(0)
129 {
131 }
132 
137 {
138 }
139 
145 {
147 
149  //Automatic determination of the analysis mode//
151  AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
152  if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
153  SetAODAnalysis();
154  } else {
155  SetESDAnalysis();
156  }
157  printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
158 
159 
164 
165  TIter next(fHistManager.GetListOfHistograms());
166  TObject* obj = 0;
167  while ((obj = next())) {
168  fOutput->Add(obj);
169  }
170  fOutput->Add(fdEdx);
171  fOutput->Add(fM20);
172  fOutput->Add(fM02);
173  fOutput->Add(fM20EovP);
174  fOutput->Add(fM02EovP);
175  fOutput->Add(fInvmassLS);
176  fOutput->Add(fInvmassULS);
177  fOutput->Add(fEMCTrketa);
178  fOutput->Add(fEMCTrkphi);
179  fOutput->Add(fHistJetEovP);
180  fOutput->Add(fHistJetEovPvPt);
181  fOutput->Add(fHistClusEovP);
184 
185 
186  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
187 }
188 
189 /*
190  * This function allocates the histograms for basic EMCal cluster QA.
191  * A set of histograms (energy, eta, phi, number of cluster) is allocated
192  * per each cluster container and per each centrality bin.
193  */
195 {
196  TString histname;
197  TString histtitle;
198  TString groupname;
199  AliClusterContainer* clusCont = 0;
200  TIter next(&fClusterCollArray);
201  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
202  groupname = clusCont->GetName();
203  fHistManager.CreateHistoGroup(groupname);
204 
205  fHistClusEovP = new TH1F("ClusEovP","Cluster EovP",130,0,1.3);
206  fHistClusEovPnonlin = new TH1F("ClusEovPnonlin","nonlin Cluster EovP",130,0,1.3);
207  fHistClusEovPHadCorr = new TH1F("ClusEovPHaddCorr","HadCorr Cluster EovP",130,0,1.3);
208  fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
209  fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
210 
211  for (Int_t cent = 0; cent < fNcentBins; cent++) {
212  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
213  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
214  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
215 
216  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
217  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
218  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
219 
220  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
221  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
222  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
223 
224  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
225  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
226  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
227 
228  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
229  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
230  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
231 
232  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
233  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
234  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
235 
236  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
237  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
238  if (fForceBeamType != kpp) {
239  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
240  }
241  else {
242  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
243  }
244  }
245  }
246 }
247 
248 /*
249  * This function allocates the histograms for basic EMCal QA.
250  * One 2D histogram with the cell energy spectra and the number of cells
251  * per event is allocated per each centrality bin.
252  */
254 {
255  TString histname;
256  TString histtitle;
257  TString groupname(fCaloCellsName);
258 
259  fHistManager.CreateHistoGroup(groupname);
260  for (Int_t cent = 0; cent < fNcentBins; cent++) {
261  histname = TString::Format("%s/histCellEnergyvsAbsId_%d", groupname.Data(), cent);
262  histtitle = TString::Format("%s;cell abs. ID;#it{E}_{cell} (GeV);counts", histname.Data());
263  fHistManager.CreateTH2(histname, histtitle, 20000, 0, 20000, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
264 
265  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
266  histtitle = TString::Format("%s;number of cells;events", histname.Data());
267  if (fForceBeamType != kpp) {
268  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
269  }
270  else {
271  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
272  }
273  }
274 }
275 
276 /*
277  * This function allocates the histograms for basic tracking QA.
278  * A set of histograms (pT, eta, phi, difference between kinematic properties
279  * at the vertex and at the EMCal surface, number of tracks) is allocated
280  * per each particle container and per each centrality bin.
281  */
283 {
284  TString histname;
285  TString histtitle;
286  TString groupname;
287  AliParticleContainer* partCont = 0;
288  TIter next(&fParticleCollArray);
289  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
290  groupname = partCont->GetName();
291  fHistManager.CreateHistoGroup(groupname);
292  for (Int_t cent = 0; cent < fNcentBins; cent++) {
293  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
294  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
295  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
296 
297  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
298  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
299  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
300 
301  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
302  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
303  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
304 
305  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
306  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
307  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
308  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
309 
310  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
311  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
312  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
313 
314  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
315  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());
316  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
317 
318  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
319  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
320  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
321  }
322 
323  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
324  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
325  if (fForceBeamType != kpp) {
326  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
327  }
328  else {
329  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
330  }
331  }
332  }
333 }
334 
335 /*
336  * This function allocates the histograms for basic jet QA.
337  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
338  * per each jet container and per each centrality bin.
339  */
341 {
342  TString histname;
343  TString histtitle;
344  TString groupname;
345  AliJetContainer* jetCont = 0;
346  TIter next(&fJetCollArray);
347  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
348  groupname = jetCont->GetName();
349  fHistManager.CreateHistoGroup(groupname);
350 
351  fHistJetEovP = new TH1F("JetEovP","HFE Jet EovP",200,0,5);
352  fHistJetEovPvPt = new TH2F("JetEovPvPt","HFE Jet EovP vs Pt",100,0,1,100,0,100);
353 
354  for (Int_t cent = 0; cent < fNcentBins; cent++) {
355  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
356  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
357  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
358 
359  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
360  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
361  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
362 
363  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
364  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
365  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
366 
367  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
368  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
369  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
370 
371  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
372  histtitle = TString::Format("%s;number of jets;events", histname.Data());
373  if (fForceBeamType != kpp) {
374  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
375  }
376  else {
377  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
378  }
379 
380  if (!jetCont->GetRhoName().IsNull()) {
381  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
382  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
383  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
384  }
385  }
386  }
387 }
388 
396 {
397  DoJetLoop();
398  DoTrackLoop();
399  DoClusterLoop();
400  DoCellLoop();
401 
402  return kTRUE;
403 }
404 
410 {
411  TString histname;
412  TString groupname;
413  AliJetContainer* jetCont = 0;
414  TIter next(&fJetCollArray);
415 
416  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
417  groupname = jetCont->GetName();
418  UInt_t count = 0;
419  Double_t rhoVal = 0;
420  rhoVal = jetCont->GetRhoVal();
421 
422  for(auto jet : jetCont->accepted()) {
423  if (!jet) continue;
424  count++;
425 
426  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
427  fHistManager.FillTH1(histname, jet->Pt());
428 
429  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
430  fHistManager.FillTH1(histname, jet->Area());
431 
432  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
433  fHistManager.FillTH1(histname, jet->Phi());
434 
435  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
436  fHistManager.FillTH1(histname, jet->Eta());
437 
438  if (jetCont->GetRhoParameter()) {
439  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
440  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
441  }
442  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
443  fHistManager.FillTH1(histname, count);
444 
445 
446 
447  Float_t ptLeading = jetCont->GetLeadingHadronPt(jet);
448  Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
449  TLorentzVector leadPart;
450  jetCont->GetLeadingHadronMomentum(leadPart, jet);
451 
452  //cout<<"JET Corr Pt: "<<corrPt<<" Pt of Leading Hadron: "<< ptLeading<<endl;
453 
454  AliParticleContainer* tracks = jetCont->GetParticleContainer();
455  if (tracks) {
456  for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
457  AliVParticle *track = jet->TrackAt(it, tracks->GetArray());
458 
459  if (track) {
460  Double_t JetTrackP = track->P();
461  Double_t JetTrackPt = track->Pt();
462 
463 
464  Double_t dphi = TVector2::Phi_0_2pi(track->Phi() - jet->Phi());
465  Double_t deta = track->Eta() - jet->Eta();
466  Double_t dist = TMath::Sqrt(deta * deta + dphi * dphi);
467 
468  //+++++----Track matching to EMCAL
469  AliVTrack *Vtrack = dynamic_cast<AliVTrack*>(track);
470  Int_t EMCalIndex = -1;
471  EMCalIndex = Vtrack->GetEMCALcluster();
472  if(EMCalIndex < 0) continue;
473 
474  AliVCluster *clustMatch = (AliVCluster*)fVevent->GetCaloCluster(EMCalIndex);
475  if(!(clustMatch && clustMatch->IsEMCAL())) continue;
476 
477  Double_t clustMatchE = clustMatch->E();
478  Double_t clustMatchNonLinE = clustMatch->GetNonLinCorrEnergy();
479  Double_t clustMatchHadCorr = clustMatch->GetHadCorrEnergy();
480 
481  //cout<<"Other Track Pt: "<< JetTrackPt << " Track Matched to Cluster Energy: "<< clustMatchE<<endl;
482  fHistJetEovP->Fill(clustMatchNonLinE / Vtrack->P());
483 
484 
485  }// Accepted Jet Track
486  }//Track loop
487  }//Track Container
488  } //jet loop
489  }//Jet Container
490 }
491 
497 {
499 
500  TString histname;
501  TString groupname;
502  AliParticleContainer* partCont = 0;
503  TIter next(&fParticleCollArray);
504  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
505  groupname = partCont->GetName();
506  UInt_t count = 0;
507  for(auto part : partCont->accepted()) {
508  if (!part) continue;
509  count++;
510 
511  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
512  fHistManager.FillTH1(histname, part->Pt());
513 
514  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
515  fHistManager.FillTH1(histname, part->Phi());
516 
517  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
518  fHistManager.FillTH1(histname, part->Eta());
519 
520  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
521  const AliVTrack* track = static_cast<const AliVTrack*>(part);
522 
523  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
524  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
525 
526  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
527  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
528 
529  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
530  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
531 
532  Double_t TrackP = track->P();
533 
535  //Photonic Electrons//
537 
538  Bool_t fFlagPhotonicElec = kFALSE;
539 
540  //----+++Identify Non-HFE
541  //SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
542 
543  if (clusCont) {
544  Int_t iCluster = track->GetEMCALcluster();
545  if (iCluster >= 0) {
546  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
547  if (cluster) {
548 
549  Double_t ClusterE = cluster->E();
550  Double_t ClusterNonLinE = cluster->GetNonLinCorrEnergy();
551  Double_t ClusterHadCorr = cluster->GetHadCorrEnergy();
552 
553  Double_t sigma = 0.0;
554 
555 
556  sigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
557  //cout<<"Sigma: " << sigma << endl;
558  fHistClusEovP->Fill(ClusterE / TrackP);
559  fHistClusEovPnonlin->Fill( ClusterNonLinE / TrackP);
560  fHistClusEovPHadCorr->Fill( ClusterHadCorr / TrackP);
561 
562  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
563  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
564  }
565  }
566  }
567  }
568  }
569 
570  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
571  fHistManager.FillTH1(histname, count);
572  }
573 }
574 
580 {
581  TString histname;
582  TString groupname;
583  AliClusterContainer* clusCont = 0;
584  TIter next(&fClusterCollArray);
585  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
586  groupname = clusCont->GetName();
587 
588  for(auto cluster : clusCont->all()) {
589  if (!cluster) continue;
590 
591  if (cluster->GetIsExotic()) {
592  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
593  fHistManager.FillTH1(histname, cluster->E());
594  }
595  }
596 
597  UInt_t count = 0;
598  for(auto cluster : clusCont->accepted()) {
599  if (!cluster) continue;
600  count++;
601 
602  AliTLorentzVector nPart;
603  cluster->GetMomentum(nPart, fVertex);
604 
605  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
606  fHistManager.FillTH1(histname, cluster->E());
607 
608  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
609  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
610 
611  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
612  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
613 
614  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
615  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
616 
617  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
618  fHistManager.FillTH1(histname, nPart.Eta());
619  }
620 
621  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
622  fHistManager.FillTH1(histname, count);
623  }
624 }
625 
631 {
632  if (!fCaloCells) return;
633 
634  TString histname;
635 
636  const Short_t ncells = fCaloCells->GetNumberOfCells();
637 
638  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
639  fHistManager.FillTH1(histname, ncells);
640 
641  histname = TString::Format("%s/histCellEnergyvsAbsId_%d", fCaloCellsName.Data(), fCentBin);
642  for (Short_t pos = 0; pos < ncells; pos++) {
643  Double_t amp = fCaloCells->GetAmplitude(pos);
644  Short_t absId = fCaloCells->GetCellNumber(pos);
645 
646  fHistManager.FillTH2(histname, absId, amp);
647  }
648 }
649 
655 {
657 }
658 
667 {
668  UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
669 
670  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
671  if (!fVevent) {
672  printf("ERROR: fVEvent not available\n");
673  return kFALSE;
674  }
675 
676  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
677  if (fESD) {
678  // printf("fESD available\n");
679  //return;
680  }
681 
682  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
683  if (fAOD) {
684  // printf("fAOD available\n");
685  //return;
686  }
687 
688  const AliVVertex *pVtx = fVevent->GetPrimaryVertex();
689 
690  fpidResponse = fInputHandler->GetPIDResponse();
691  if (!fpidResponse) {
692  printf("ERROR: fpidResponse not available\n");
693  return kFALSE;
694  }
695 
696  Int_t ntracks;
697  ntracks = fVevent->GetNumberOfTracks();
698  //printf("There are %d tracks in this event\n",ntracks);
699 
700  Double_t Zvertex = -100, Xvertex = -100, Yvertex = -100;
701  Double_t NcontV = pVtx->GetNContributors();
702 
703  // if(NcontV<2)return kTRUE; //Events with 2 Tracks
704 
706  //EMcal cluster info//
709 
710 
711  Int_t numberofvertices = 100;
712  if(fAOD) numberofvertices = fAOD->GetNumberOfVertices();
713  Double_t listofmotherkink[numberofvertices];
714  Int_t numberofmotherkink = 0;
715  if(IsAODanalysis())
716  {
717  for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
718  AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
719  if(!aodvertex) continue;
720  if(aodvertex->GetType()==AliAODVertex::kKink) {
721  AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
722  if(!mother) continue;
723  Int_t idmother = mother->GetID();
724  listofmotherkink[numberofmotherkink] = idmother;
725  numberofmotherkink++;
726  }
727  }
728  } //+++
729 
730  fEventCounter++;
731 
732  //cout<<"Event: "<<fEventCounter<<" Number of Vertices: "<<numberofvertices<<" Number of Tracks: "<<ntracks<< endl;
733 
734  return kTRUE;
735 }
736 
737 //________________________________________________________________________
739 {
741  //EMCAL cluster information//
743 
744  Int_t Nclust = -999;
745  TVector3 clustpos;
746  Float_t emcx[3]; // cluster pos
747  Double_t clustE=-999, emcphi = -999, emceta=-999;
748  Nclust = fVevent->GetNumberOfCaloClusters();
749  for(Int_t icl=0; icl<Nclust; icl++)
750  {
751  AliVCluster *clust = 0x0;
752  clust = fVevent->GetCaloCluster(icl);
753  if(!clust) printf("ERROR: Could not receive cluster matched calibrated from track %d\n", icl);
754 
755  if(clust && clust->IsEMCAL())
756  {
757  clustE = clust->E();
758  clust->GetPosition(emcx);
759  clustpos.SetXYZ(emcx[0],emcx[1],emcx[2]);
760  emcphi = clustpos.Phi();
761  emceta = clustpos.Eta();
762  //fHistClustE->Fill(clustE);
763  //fEMCClsEtaPhi->Fill(emceta,emcphi);
764  }
765  }
766 }
767 //________________________________________________________________________
768 
769 //___________________________________________________________________________
770 
771 void AliAnalysisTaskEmcalJetHF::SelectPhotonicElectron(Int_t itrack, AliVTrack *track, Bool_t &fFlagPhotonicElec)
772 {
774  //Photonic electron selection//
776 
777  Bool_t flagPhotonicElec = kFALSE;
778  Double_t ptAsso=-999., nsigma=-999.0;
779  Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
780 
781  for(Int_t jTracks = 0; jTracks<fVevent->GetNumberOfTracks(); jTracks++){
782  if(jTracks==itrack) continue;
783 
784  AliVParticle* VtrackAsso = fVevent->GetTrack(jTracks);
785  if (!VtrackAsso) {
786  printf("ERROR: Could not receive track %d\n", jTracks);
787  continue;
788  }
789 
790  AliVTrack *trackAsso = dynamic_cast<AliVTrack*>(VtrackAsso);
791  if(!trackAsso) continue;
792 
793  if(IsAODanalysis()) {
794  AliAODTrack *atrackAsso = dynamic_cast<AliAODTrack*>(VtrackAsso);
795  if(!atrackAsso) continue;
796  if(!atrackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
797  if(atrackAsso->GetTPCNcls() < 70) continue;
798  if(!(atrackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
799  if(!(atrackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
800  }
801 
802  nsigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
803  ptAsso = trackAsso->Pt();
804  Int_t chargeAsso = trackAsso->Charge();
805  Int_t charge = track->Charge();
806 
807  if(ptAsso <0.2) continue;
808  if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
809  if(nsigma < -3 || nsigma > 3) continue;
810 
811  Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
812  if(charge>0) fPDGe1 = -11;
813  if(chargeAsso>0) fPDGe2 = -11;
814 
815  if(charge == chargeAsso) fFlagLS = kTRUE;
816  if(charge != chargeAsso) fFlagULS = kTRUE;
817 
818  AliKFParticle::SetField(fVevent->GetMagneticField());
819 
820  AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
821  AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
822  AliKFParticle recg(ge1, ge2);
823 
824  if(recg.GetNDF()<1) continue;
825  Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
826  if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
827 
828  Double_t mass=-999., width = -999;
829  Int_t MassCorrect;
830  MassCorrect = recg.GetMass(mass,width);
831 
832  if(fFlagLS && track->Pt()>1) fInvmassLS->Fill(mass);
833  if(fFlagULS && track->Pt()>1) fInvmassULS->Fill(mass);
834 
835  if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
836  flagPhotonicElec = kTRUE;
837  }
838  }
839  fFlagPhotonicElec = flagPhotonicElec;
840 }
841 
842 
847 {
848  cout<<"#######################"<<endl;
849  cout<<"#### Task Finished ####"<<endl;
850  cout<<"#######################"<<endl;
851 }
Int_t charge
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
const TString & GetRhoName() const
AliAODEvent * fAOD
ESD object.
Double_t mass
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
Int_t fCentBin
!event centrality bin
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
THashList * CreateHistoGroup(const char *groupname, const char *parent="/")
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="")
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
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
Implementation of a HFE jet analysis task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
TObjArray fJetCollArray
jet collection array
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
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 AliJetIterableContainer accepted() const
Container structure for EMCAL clusters.
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
Declaration of class AliAnalysisTaskEmcalJetHF.