AliPhysics  aaf9c62 (aaf9c62)
AliAnalysisTaskChargedParticlesMCTriggerMimic.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2015, 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 #include <array>
16 #include <bitset>
17 #include <iomanip>
18 #include <iostream>
19 #include <memory>
20 #include <string>
21 
22 #include <TClonesArray.h>
23 #include <TFile.h>
24 #include <THistManager.h>
25 #include <TLinearBinning.h>
26 #include <TMath.h>
27 #include <TPDGCode.h>
28 #include <TTree.h>
29 
30 #include "AliAnalysisUtils.h"
31 #include "AliAODInputHandler.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODTrack.h"
35 #include "AliClusterContainer.h"
37 #include "AliEMCALGeometry.h"
38 #include "AliEMCALRecoUtils.h"
39 #include "AliEMCALTriggerPatchInfo.h"
41 #include "AliESDtrack.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "AliInputEventHandler.h"
44 #include "AliLog.h"
45 #include "AliVCluster.h"
46 #include "AliVEvent.h"
47 #include "AliVParticle.h"
48 
50 
54 
56 
60 AliAnalysisTaskChargedParticlesMCTriggerMimic::AliAnalysisTaskChargedParticlesMCTriggerMimic():
62  fTrackCuts(NULL),
63  fHistos(NULL),
64  fWeightHandler(NULL),
65  fYshift(0.465),
66  fEtaSign(1),
67  fEtaLabCut(-0.5, 0.5),
68  fEtaCmsCut(-0.13, 0.13),
69  fPhiCut(0, TMath::TwoPi()),
70  fPatchType(kUndef),
71  fEnergyThreshold(0.),
72  fObservables(),
73  fNameClusters()
74 {
75  SetCaloTriggerPatchInfoName("EmcalTriggers");
76  SetNeedEmcalGeom(true);
77 }
78 
84  AliAnalysisTaskEmcal(name, true),
88  fYshift(0.465),
89  fEtaSign(1),
90  fEtaLabCut(-0.6, 0.6),
91  fEtaCmsCut(-0.13, 0.13),
93  fEnergyThreshold(0.),
94  fObservables(),
96 {
97  SetCaloTriggerPatchInfoName("EmcalTriggers");
98  SetNeedEmcalGeom(true);
100 }
101 
106  //if(fTrackCuts) delete fTrackCuts;
107  if(fHistos) delete fHistos;
108 }
109 
112 
113  if(fIsPythia){
114  AliDebugStream(1) << GetName() << ": Running on PYTHIA Hard production" << std::endl;
115  } else {
116  AliDebugStream(1) << GetName() << ": Not running on PYTHIA Hard production" << std::endl;
117  }
118 
119  if(!fAliAnalysisUtils) fAliAnalysisUtils = new AliAnalysisUtils;
120  if(!fTrackCuts) InitializeTrackCuts("standard",fInputHandler->IsA() == AliAODInputHandler::Class());
121 
122  PtBinning newbinning;
123  TLinearBinning smbinning(21, -0.5, 20.5), etabinning(100, -0.7, 0.7);
124  fHistos = new THistManager("Ref");
125 
126  // The first two histograms should lead to a re-implementation inside AliAnalysisTaskEMCAL:
127  // Users should be able to set event weights.
128  fHistos->CreateTH1("hUserEventCount", "User event counter", 1, 0.5, 1.5);
129  fHistos->CreateTH1("hUserVertexZ", "User vertex distribution after z-cut", 100, -10, 10);
130  fHistos->CreateTH1("hUserPtHard", "User pt-hard distribution", 1000, 0., 300.);
131 
132  // Histograms for observable tracks
133  if(HasObservable(kTracks)){
134  const std::array<std::string, 2> kInputs = {"True", "Accept"};
135  const std::array<std::string, 6> kSpecies = {"El", "Mu", "Pi", "Ka", "Pr", "Ot"};
136  const std::array<double, 5> kPtCuts = {1., 2., 5., 10., 20.};
137  for(const auto &input : kInputs){
138  AliDebugStream(1) << GetName() << ": Creating histograms for case " << input << std::endl;
139  fHistos->CreateTH1(Form("hTrackPtEtaAll%s", input.c_str()), Form("Charged particle p_{t} distribution all #eta %s", input.c_str()), newbinning, "s");
140  fHistos->CreateTH1(Form("hTrackPtEtaCent%s", input.c_str()), Form("Charged particle p_{t} distribution central #eta %s", input.c_str()), newbinning, "s");
141  fHistos->CreateTH1(Form("hTrackPtEMCALEtaAll%s", input.c_str()), Form("Charged particle in EMCAL p_{t} distribution all #eta trigger %s", input.c_str()), newbinning);
142  fHistos->CreateTH1(Form("hTrackPtEMCALEtaCent%s", input.c_str()), Form("Charged particle in EMCAL p_{t} distribution central eta trigger %s", input.c_str()), newbinning);
143  for(const auto &piditer : kSpecies){
144  fHistos->CreateTH1(Form("hTrackPtEtaAll%s%s", piditer.c_str(), input.c_str()), Form("Charged %s p_{t} distribution all #eta %s", piditer.c_str(), input.c_str()), newbinning);
145  fHistos->CreateTH1(Form("hTrackPtEtaCent%s%s", piditer.c_str(), input.c_str()), Form("Charged %s p_{t} distribution central #eta %s", piditer.c_str(), input.c_str()), newbinning);
146  fHistos->CreateTH1(Form("hTrackPtEMCALEtaAll%s%s", piditer.c_str(), input.c_str()), Form("Charged %s in EMCAL p_{t} distribution all #eta %s", piditer.c_str(), input.c_str()), newbinning);
147  fHistos->CreateTH1(Form("hTrackPtEMCALEtaCent%s%s", piditer.c_str(), input.c_str()), Form("Charged %s in EMCAL p_{t} distribution central #eta %s", piditer.c_str(), input.c_str()), newbinning);
148  }
149  for(const auto &ptcut : kPtCuts){
151  Form("hTrackEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
152  Form("#eta_{lab} distribution without #eta-cut for tracks with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
153  100,
154  -1.,
155  1.
156  );
158  Form("hTrackEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
159  Form("#eta_{lab} distribution with #eta-cut for tracks with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
160  100,
161  -1.,
162  1.
163  );
165  Form("hTrackEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
166  Form("#eta_{cent} distribution without #eta-cut for tracks with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
167  160,
168  -1.3,
169  1.3
170  );
172  Form("hTrackEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
173  Form("#eta_{cent} distribution with #eta-cut for tracks with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
174  160,
175  -1.3,
176  1.3
177  );
179  Form("hTrackEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
180  Form("#eta_{lab} distribution without #eta-cut for tracks in EMCAL with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
181  100,
182  -1.,
183  1.
184  );
186  Form("hTrackEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
187  Form("#eta_{lab} distribution with #eta-cut for tracks in EMCAL with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
188  100,
189  -1.,
190  1.
191  );
193  Form("hTrackEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
194  Form("#eta_{cent} distribution without #eta-cut for tracks in EMCAL with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
195  160,
196  -1.3,
197  1.3
198  );
200  Form("hTrackEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
201  Form("Eta (cent) distribution with #eta-cut for tracks in EMCAL with p_{t} above %.1f GeV/c %s", ptcut, input.c_str()),
202  160,
203  -1.3,
204  1.3
205  );
207  Form("hTrackPhiDistAllPt%d%s", static_cast<Int_t>(ptcut), input.c_str()),
208  Form("#phi distribution of particles with p_{t} above %.1f GeV/c trigger %s", ptcut, input.c_str()),
209  300,
210  0.,
211  2*TMath::Pi()
212  );
213  }
214  }
215  }
216 
217  // Histograms for observable clusters
218  std::array<Double_t, 5> kEnCuts = {1., 2., 5., 10., 20.};
220  Int_t sectorsWithEMCAL[10] = {4, 5, 6, 7, 8, 9, 13, 14, 15, 16};
221  fHistos->CreateTH1("hClusterEnergy", "Cluster energy", newbinning);
222  fHistos->CreateTH1("hClusterET", "Cluster transverse energy", newbinning);
223  fHistos->CreateTH2("hClusterEnergySM", "Cluster energy versus supermodule", smbinning, newbinning);
224  fHistos->CreateTH2("hClusterETSM", "Cluster transverse energy versus supermodule", smbinning, newbinning);
225  fHistos->CreateTH2("hClusterEtaEnergy", "Cluster energy vs. eta", etabinning, newbinning);
226  fHistos->CreateTH2("hClusterEtaET", "Cluster transverse energy vs. eta", etabinning, newbinning);
227  for(int ism = 0; ism < 20; ism++){
228  fHistos->CreateTH2(Form("hClusterEtaEnergySM%d", ism), Form("Cluster energy vs. eta in Supermodule %d", ism), etabinning, newbinning);
229  fHistos->CreateTH2(Form("hClusterEtaETSM%d", ism), Form("Cluster transverse energy vs. eta in Supermodule %d", ism), etabinning, newbinning);
230  }
231  for(int isec = 0; isec < 10; isec++){
232  fHistos->CreateTH2(Form("hClusterEtaEnergySec%d", sectorsWithEMCAL[isec]), Form("Cluster energy vs.eta in tracking sector %d", sectorsWithEMCAL[isec]), etabinning, newbinning);
233  fHistos->CreateTH2(Form("hClusterEtaETSec%d", sectorsWithEMCAL[isec]), Form("Cluster transverse energy vs.eta in tracking sector %d", sectorsWithEMCAL[isec]), etabinning, newbinning);
234  }
235  for(auto ien : kEnCuts){
236  fHistos->CreateTH2(Form("hClusterEtaPhi%dG", static_cast<int>(ien)), Form("cluster #eta-#phi map for clusters with energy larger than %f GeV/c", ien), 100, -0.7, 0.7, 200, 0, 2*TMath::Pi());
237  }
238  }
239 
240  // Histograms for Observable EGA or EJE patches
242  const int kNPatchTypes = 2;
243  const std::array<Observable_t, kNPatchTypes> kPatchObservables = {kEGAPatches, kEJEPatches};
244  const std::array<TString, kNPatchTypes> kPatchTypes = {"EGA","EJE"};
245  for(int ipatch = 0; ipatch < kNPatchTypes; ipatch++){
246  if(!HasObservable(kPatchObservables[ipatch])) continue;
247  fHistos->CreateTH1(Form("h%sPatchEnergy", kPatchTypes[ipatch].Data()), Form("%s-patch energy", kPatchTypes[ipatch].Data()), newbinning);
248  fHistos->CreateTH1(Form("h%sPatchET", kPatchTypes[ipatch].Data()), Form("%s-patch transverse energy", kPatchTypes[ipatch].Data()), newbinning);
249  fHistos->CreateTH2(Form("h%sPatchEnergyEta", kPatchTypes[ipatch].Data()), Form("%s-patch energy", kPatchTypes[ipatch].Data()), newbinning, etabinning);
250  fHistos->CreateTH2(Form("h%sPatchETEta", kPatchTypes[ipatch].Data()), Form("%s-patch transverse energy", kPatchTypes[ipatch].Data()), newbinning, etabinning);
251  for(auto enc : kEnCuts){
252  fHistos->CreateTH2(Form("h%sPatchEtaPhi%dG", kPatchTypes[ipatch].Data(), static_cast<int>(enc)), Form("%s-patch #eta-#phi map for patches with energy larger than %f GeV/c", kPatchTypes[ipatch].Data(), enc), 100, -0.7, 0.7, 200, 0, TMath::TwoPi());
253  fHistos->CreateTH2(Form("h%sPatchColRow%dG", kPatchTypes[ipatch].Data(), static_cast<int>(enc)), Form("%s-patch col-row map for patches with energy larger than %f GeV/c", kPatchTypes[ipatch].Data(), enc), 48, -0.5, 47.5, 104, -0.5, 103.5);
254  }
255  }
256  }
257 
258  for(auto hist : *(fHistos->GetListOfHistograms())) fOutput->Add(hist);
259  PostData(1, fOutput);
260 
261  AliDebugStream(1) << GetName() << ": Output objects initialized" << std::endl;
262 }
263 
270  AliDebugStream(2) << GetName() << ": Using custom event selection method" << std::endl;
271  if(!fTriggerPatchInfo){
272  AliErrorStream() << GetName() << ": Trigger patch container not found but required" << std::endl;
273  return false;
274  }
275  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
276  AliDebugStream(3) << GetName() << "Event is an INT7 event" << std::endl;
277 
278  // MC outlier cut
279  if(fIsPythia){
280  if(!CheckMCOutliers()){
281  AliDebugStream(3) << GetName() << ": Event identified as outlier" << std::endl;
282  return false;
283  } else {
284  AliDebugStream(3) << GetName() << ": Not an outlier event" << std::endl;
285  }
286  }
287 
288  // Generall event quality cuts
289  // The vertex cut also contains cuts on the number
290  // of contributors and the position in z
291  AliDebugStream(3) << GetName() << ": Applying vertex selection" << std::endl;
292  if(fAliAnalysisUtils){
293  if(!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) return false;
294  if(fAliAnalysisUtils->IsPileUpEvent(InputEvent())) return false;
295  AliDebugStream(3) << GetName() << ": Vertex selection passed" << std::endl;
296  }
297 
298  AliDebugStream(3) << GetName() << ": Applying EMCAL trigger selection" << std::endl;
299  if(fPatchType != kUndef){
301  AliDebugStream(3) << GetName() << ": Failed trigger selection" << std::endl;
302  return false;
303  }
304  }
305 
306  AliDebugStream(2) << GetName() << "Event selected" << std::endl;
307  return true;
308 }
309 
315  AliDebugStream(1) << GetName() << ": Inspecting event" << std::endl;
316  Double_t weight = 1.;
317  if(fIsPythia){
318  if(!fPythiaHeader){
319  AliErrorStream() << GetName() << ": PYTHIA event header not found" << std::endl;
320  } else {
322  }
323  }
324 
325  fHistos->FillTH1("hUserEventCount", 1, weight);
326  const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
327  fHistos->FillTH1("hUserVertexZ", vtx->GetZ(), weight);
328  fHistos->FillTH1("hUserPtHard", fPtHard);
329 
330  if(HasObservable(kTracks)){
331  AliDebugStream(3) << GetName() << ": eta-lab cut: " << fEtaLabCut << ", eta-cms cut: " << fEtaCmsCut << std::endl;
332 
333  AliVParticle *truepart(nullptr);
334  Bool_t isEMCAL(kFALSE);
335  if(MCEvent()){
336  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
337  truepart = MCEvent()->GetTrack(ipart);
338 
339  // Select only particles within ALICE acceptance
340  if(!fEtaLabCut.IsInRange(truepart->Eta())) continue;
341  if(!fPhiCut.IsInRange(truepart->Phi())) continue;
342  if(TMath::Abs(truepart->Pt()) < 0.1) continue;
343  if(!truepart->Charge()) continue;
344 
345  if(!IsPhysicalPrimary(truepart, fMCEvent)) continue;
346  AliAODMCParticle *aodmc = static_cast<AliAODMCParticle *>(truepart);
347  isEMCAL = (truepart->Phi() > 1.5 && truepart->Phi() < 3.1) ? kTRUE : kFALSE;
348 
349  // Calculate eta in cms frame according
350  // EPJC74 (2014) 3054:
351  // eta_cms = - eta_lab - |yshift|
352  Double_t etacent = -1. * truepart->Eta() - TMath::Abs(fYshift);
353  etacent *= fEtaSign;
354 
355  Bool_t etacentcut = fEtaCmsCut.IsInRange(etacent);
356 
357  // Get PID
358  TString pid = "";
359  switch(TMath::Abs(truepart->PdgCode())){
360  case kPiPlus: pid = "Pi"; break;
361  case kMuonMinus: pid = "Mu"; break;
362  case kElectron: pid = "El"; break;
363  case kKPlus: pid = "Ka"; break;
364  case kProton: pid = "Pr"; break;
365  default: pid = "Ot"; break;
366  };
367 
368  // Particle selected (do not filter TRD sectors for MC truth)
369  FillTrackHistos("True", weight, TMath::Abs(truepart->Pt()), truepart->Eta() * fEtaSign, etacent, truepart->Phi(), etacentcut, isEMCAL, pid);
370  }
371  }
372 
373  // Loop over tracks, fill select particles
374  // Histograms
375  // - Full eta (-0.8, 0.8), new binning
376  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
377  // - Central eta (-0.8, -0.2), new binning,
378  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
379  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
380  AliVTrack *checktrack(NULL);
381  AliVParticle *assocMC(NULL);
382  double ptparticle(-1.), etaparticle(-100.), etaEMCAL(0.), phiEMCAL(0.);
383  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
384  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
385  if(!checktrack) continue;
386  TString assocpid = "Ot";
387  // Find associated particle
388  if(MCEvent()){
389  assocMC = MCEvent()->GetTrack(TMath::Abs(checktrack->GetLabel()));
390  if(!assocMC) continue; // Fake track
391  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
392  // Get PID
393  switch(TMath::Abs(assocMC->PdgCode())){
394  case kPiPlus: assocpid = "Pi"; break;
395  case kMuonMinus: assocpid = "Mu"; break;
396  case kElectron: assocpid = "El"; break;
397  case kKPlus: assocpid = "Ka"; break;
398  case kProton: assocpid = "Pr"; break;
399  default: assocpid = "Ot"; break;
400  };
401  }
402 
403  // Select only particles within ALICE acceptance
404  if(!fEtaLabCut.IsInRange(checktrack->Eta())) continue;
405  if(!fPhiCut.IsInRange(checktrack->Phi())) continue;
406  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
407  if(checktrack->IsA() == AliESDtrack::Class()){
408  AliESDtrack copytrack(*(static_cast<AliESDtrack *>(checktrack)));
409  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
410  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
411  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
412  } else {
413  AliAODTrack copytrack(*(static_cast<AliAODTrack *>(checktrack)));
414  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
415  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
416  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
417  }
418  Int_t supermoduleID = -1;
419  isEMCAL = fGeom->SuperModuleNumberFromEtaPhi(etaEMCAL, phiEMCAL, supermoduleID);
420  // Exclude supermodules 10 and 11 as they did not participate in the trigger
421  isEMCAL = isEMCAL && supermoduleID < 10;
422 
423  if(!fTrackCuts->IsTrackAccepted(checktrack)) continue;
424 
425  // prefer true pt and eta, however in case of running on data take measured values
426  ptparticle = assocMC ? TMath::Abs(assocMC->Pt()) : TMath::Abs(checktrack->Pt());
427  etaparticle = assocMC ? assocMC->Eta() : checktrack->Eta();
428 
429  // Calculate eta in cms frame according
430  // EPJC74 (2014) 3054:
431  // eta_cms = - eta_lab - |yshift|
432  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
433  etacent *= fEtaSign;
434 
435  Bool_t etacentcut = fEtaCmsCut.IsInRange(etacent);
436 
437  FillTrackHistos("Accept", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, assocpid);
438  }
439  }
440 
442  Double_t vertexpos[3];
443  fInputEvent->GetPrimaryVertex()->GetXYZ(vertexpos);
444 
445  Double_t energy, et, eta, phi;
446  AliClusterContainer *clustercont = this->GetClusterContainer(fNameClusters.Data());
447  if(clustercont){
448  for(const auto &clust : clustercont->all()){
449  if(!clust->IsEMCAL()) continue;
450  if(clust->GetIsExotic()) continue;
451 
452  TLorentzVector posvec;
453  energy = clust->GetNonLinCorrEnergy();
454  et = posvec.Et();
455  clust->GetMomentum(posvec, vertexpos);
456  eta = posvec.Eta();
457  phi = posvec.Phi();
458  FillClusterHistos(weight, energy, et, eta, phi);
459  }
460  }
461  }
462 
464  AliEMCALTriggerPatchInfo *recpatch(nullptr);
465  TString patchname;
466  for(TIter patchiter = TIter(fTriggerPatchInfo).Begin(); patchiter != TIter::End(); ++patchiter){
467  recpatch = static_cast<AliEMCALTriggerPatchInfo *>(*patchiter);
468  if(!recpatch->IsOfflineSimple()) continue;
469  if(recpatch->IsJetHighSimple() && HasObservable(kEJEPatches)){
470  patchname = "EJE";
471  } else if(recpatch->IsGammaHighSimple() && HasObservable(kEGAPatches)){
472  patchname = "EGA";
473  } else continue;
474 
475  FillPatchHistos(patchname.Data(), weight, recpatch->GetPatchE(), recpatch->GetPatchET(), recpatch->GetEtaGeo(), recpatch->GetPhiGeo(), recpatch->GetColStart(), recpatch->GetRowStart());
476  }
477  }
478 
479  return true;
480 }
481 
494  const char *eventclass,
495  Double_t weight,
496  Double_t pt,
497  Double_t etalab,
498  Double_t etacent,
499  Double_t phi,
500  Bool_t etacut,
501  Bool_t inEmcal,
502  const char *pid
503  )
504 {
505  fHistos->FillTH1(Form("hTrackPtEtaAll%s", eventclass), TMath::Abs(pt), weight);
506  fHistos->FillTH1(Form("hTrackPtEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
507  if(inEmcal){
508  fHistos->FillTH1(Form("hTrackPtEMCALEtaAll%s", eventclass), TMath::Abs(pt), weight);
509  fHistos->FillTH1(Form("hTrackPtEMCALEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
510  }
511 
512  const std::array<int, 5> kPtMin = {1,2,5,10,20}; // for eta distributions
513  for(const auto &ptmin : kPtMin){
514  if(TMath::Abs(pt) > static_cast<double>(ptmin)){
515  fHistos->FillTH1(Form("hTrackPhiDistAllPt%d%s", ptmin, eventclass), phi, weight);
516  fHistos->FillTH1(Form("hTrackEtaLabDistAllPt%d%s", ptmin, eventclass), etalab, weight);
517  fHistos->FillTH1(Form("hTrackEtaCentDistAllPt%d%s", ptmin, eventclass), etacent, weight);
518  if(inEmcal){
519  fHistos->FillTH1(Form("hTrackEtaLabDistAllEMCALPt%d%s", ptmin, eventclass), etalab, weight);
520  fHistos->FillTH1(Form("hTrackEtaCentDistAllEMCALPt%d%s", ptmin, eventclass), etacent, weight);
521  }
522  }
523  }
524 
525  if(etacut){
526  fHistos->FillTH1(Form("hTrackPtEtaCent%s", eventclass), TMath::Abs(pt), weight);
527  fHistos->FillTH1(Form("hTrackPtEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
528  if(inEmcal){
529  fHistos->FillTH1(Form("hTrackPtEMCALEtaCent%s", eventclass), TMath::Abs(pt), weight);
530  fHistos->FillTH1(Form("hTrackPtEMCALEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
531  }
532  for(const auto &ptmin : kPtMin){
533  if(TMath::Abs(pt) > static_cast<double>(ptmin)){
534  fHistos->FillTH1(Form("hTrackEtaLabDistCutPt%d%s", ptmin, eventclass), etalab, weight);
535  fHistos->FillTH1(Form("hTrackEtaCentDistCutPt%d%s", ptmin, eventclass), etacent, weight);
536  if(inEmcal){
537  fHistos->FillTH1(Form("hTrackEtaLabDistCutEMCALPt%d%s", ptmin, eventclass), etalab, weight);
538  fHistos->FillTH1(Form("hTrackEtaCentDistCutEMCALPt%d%s", ptmin, eventclass), etacent, weight);
539  }
540  }
541  }
542  }
543 }
544 
554  double weight,
555  double energy,
556  double transverseenergy,
557  double eta,
558  double phi
559  )
560 {
561  Int_t supermoduleID = -1, sector = -1;
562  fGeom->SuperModuleNumberFromEtaPhi(eta, phi, supermoduleID);
563 
564  fHistos->FillTH1("hClusterEnergy", energy, weight);
565  fHistos->FillTH1("hClusterET", transverseenergy, weight);
566  fHistos->FillTH2("hClusterEtaEnergy", eta, energy, weight);
567  fHistos->FillTH2("hClusterEtaET", eta, transverseenergy, weight);
568  if(supermoduleID >= 0){
569  fHistos->FillTH2("hClusterEnergySM", supermoduleID, energy, weight);
570  fHistos->FillTH2("hClusterETSM", supermoduleID, transverseenergy, weight);
571  fHistos->FillTH2(Form("hClusterEtaEnergySM%d", supermoduleID), eta, energy, weight);
572  fHistos->FillTH2(Form("hClusterEtaETSM%d", supermoduleID), eta, transverseenergy, weight);
573  if(supermoduleID < 12)
574  sector = 4 + int(supermoduleID/2); // EMCAL
575  else
576  sector = 13 + int((supermoduleID-12)/2); // DCAL
577  fHistos->FillTH2(Form("hClusterEtaEnergySec%d", sector), eta, energy, weight);
578  fHistos->FillTH2(Form("hClusterEtaETSec%d", sector), eta, transverseenergy, weight);
579  }
580  std::array<Double_t, 5> encuts = {1., 2., 5., 10., 20.};
581  for(auto e : encuts){
582  if(energy > e){
583  fHistos->FillTH2(Form("hClusterEtaPhi%dG", static_cast<int>(e)), eta, phi, weight);
584  }
585  }
586 }
587 
595 void AliAnalysisTaskChargedParticlesMCTriggerMimic::FillPatchHistos(const char *patchname, double weight, double energy, double transverseenergy, double eta, double phi, int col, int row){
596  fHistos->FillTH1(Form("h%sPatchEnergy", patchname), energy, weight);
597  fHistos->FillTH1(Form("h%sPatchET", patchname), transverseenergy, weight);
598  fHistos->FillTH2(Form("h%sPatchEnergyEta", patchname), energy, eta, weight);
599  fHistos->FillTH2(Form("h%sPatchETEta", patchname), transverseenergy, eta, weight);
600  const std::array<Double_t, 5> kEnCuts = {1., 2., 5., 10., 20.};
601  for(auto e : kEnCuts){
602  if(energy > e){
603  fHistos->FillTH2(Form("h%sPatchEtaPhi%dG", patchname, static_cast<int>(e)), eta, phi, weight);
604  fHistos->FillTH2(Form("h%sPatchColRow%dG", patchname, static_cast<int>(e)), col, row, weight);
605  }
606  }
607 }
608 
609 
618 Bool_t AliAnalysisTaskChargedParticlesMCTriggerMimic::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) {
619  Bool_t physprim = false;
620  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
621  if(aodmc){
622  physprim = aodmc->IsPhysicalPrimary();
623  } else {
624  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
625  }
626  return physprim;
627 }
628 
636 }
637 
646  bool selected = false;
647  AliEMCALTriggerPatchInfo *patch(nullptr);
648  AliDebugStream(2) << GetName() << ": Selecting EMCAL triggered event type (" << (fPatchType == kEMCEGA ? "EGA" : "EJE") << ") using patch energy above threshold" << std::endl;
649  AliDebugStream(2) << GetName() << ": Energy threshold " << fEnergyThreshold << " GeV" << std::endl;
650  AliDebugStream(2) << GetName() << ": Number of reconstructed patches " << triggerpatches->GetEntries() << std::endl;
651  for(TIter patchiter = TIter(triggerpatches).Begin(); patchiter != TIter::End(); ++patchiter){
652  patch = static_cast<AliEMCALTriggerPatchInfo *>(*patchiter);
653  AliDebugStream(4) << GetName() << ": Next patch" << std::endl;
654  if(!patch->IsOfflineSimple()) continue;
655  AliDebugStream(4) << GetName() << "Patch is an offline simple patch" << std::endl;
656  AliDebugStream(4) << GetName() << ": Trigger bits: " << std::bitset<32>(patch->GetTriggerBits()) << std::endl;
657  AliDebugStream(4) << GetName() << ": J1(" << patch->IsJetHighSimple() << "), J2(" << patch->IsJetLowSimple()
658  << "), G1(" << patch->IsGammaHighSimple() << ") G2(" << patch->IsGammaLowSimple() << ")" << std::endl;
659  if(fPatchType == kEMCEJE){
660  if(!patch->IsJetHighSimple()) continue;
661  AliDebugStream(3) << GetName() << ": Patch is jet high simple" << std::endl;
662  } else if(fPatchType == kEMCEGA){
663  if(!patch->IsGammaHighSimple()) continue;
664  AliDebugStream(4) << GetName() << ": Patch is gamma high simple" << std::endl;
665  }
666  AliDebugStream(3) << GetName() << ": Found trigger patch of matching type, now cutting on energy ...." << std::endl;
667  if(patch->GetPatchE() > fEnergyThreshold){
668  // firing trigger patch found
669  AliDebugStream(2) << GetName() << ": Firing trigger patch found at energy " << std::setprecision(1) << patch->GetPatchE() << std::endl;
670  selected = true;
671  break;
672  }
673  }
674  return selected;
675 }
676 
682 {
683  this->SetMinimum(0.);
684  this->AddStep(1, 0.05);
685  this->AddStep(2, 0.1);
686  this->AddStep(4, 0.2);
687  this->AddStep(7, 0.5);
688  this->AddStep(16, 1);
689  this->AddStep(36, 2);
690  this->AddStep(40, 4);
691  this->AddStep(50, 5);
692  this->AddStep(100, 10);
693  this->AddStep(200, 20);
694 }
695 
696 
697 } /* namespace EMCalTriggerPtAnalysis */
Bool_t fIsPythia
trigger, if it is a PYTHIA production
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
void FillTrackHistos(const char *eventclass, Double_t weight, Double_t pt, Double_t eta, Double_t etacent, Double_t phi, Bool_t etacut, Bool_t inEmcal, const char *pid)
Base task in the EMCAL framework.
energy
Definition: HFPtSpectrum.C:44
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
void AddStep(Double_t max, Double_t binwidth)
void FillPatchHistos(const char *patchname, double weight, double energy, double transverseenergy, double eta, double phi, int col, int row)
void SetCaloTriggerPatchInfoName(const char *n)
static AliEmcalTrackSelection * TrackCutsFactory(TString name, Bool_t isAOD)
Fully-configure EMCAL track selection independent of the data type.
const AliClusterIterableContainer all() const
Bool_t CheckMCOutliers()
Filter the mc tails in pt-hard distributions.
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
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
AliEMCALGeometry * fGeom
!emcal geometry
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliAnalysisUtils * fAliAnalysisUtils
!vertex selection (optional)
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Helper class creating user defined custom binning.
const Double_t ptmin
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Float_t fPtHard
!event -hard
double GetEventWeight(const AliMCEvent *const event) const
virtual PWG::EMCAL::AliEmcalTrackSelResultPtr IsTrackAccepted(AliVTrack *const trk)=0
Interface for track selection code.
AliEmcalList * fOutput
!output list
void FillClusterHistos(double weight, double energy, double transversenergy, double eta, double phi)
Bool_t IsPhysicalPrimary(const AliVParticle *const part, AliMCEvent *const mcevent)
Analysis of high- tracks in triggered events.
void SetMakeGeneralHistograms(Bool_t g)
TClonesArray * fTriggerPatchInfo
!trigger patch info array
void SetNeedEmcalGeom(Bool_t n)
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
Container class for histograms.
Definition: THistManager.h:99
Double_t fEtaSign
Sign of the eta distribution (swaps when beam directions swap): p-Pb: +1, Pb-p: -1.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Container structure for EMCAL clusters.
void SetMinimum(Double_t min)