AliPhysics  1a228f7 (1a228f7)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskChargedParticlesRefMC.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 <iostream>
16 #include <memory>
17 
18 #include <TArrayD.h>
19 #include <TChain.h>
20 #include <TClonesArray.h>
21 #include <TFile.h>
22 #include <THashList.h>
23 #include <TH1.h>
24 #include <THistManager.h>
25 #include <TKey.h>
26 #include <TList.h>
27 #include <TPDGCode.h>
28 #include <TProfile.h>
29 #include <TMath.h>
30 #include <TString.h>
31 #include <TSystem.h>
32 #include <TTree.h>
33 
34 #include "AliAnalysisManager.h"
35 #include "AliAnalysisUtils.h"
36 #include "AliAODMCHeader.h"
37 #include "AliAODInputHandler.h"
38 #include "AliAODMCParticle.h"
39 #include "AliAODTrack.h"
40 #include "AliEmcalAnalysisFactory.h"
41 #include "AliEmcalTrackSelection.h"
42 #include "AliEmcalTriggerOfflineSelection.h"
43 #include "AliEMCALTriggerPatchInfo.h"
44 #include "AliEMCALGeometry.h"
45 #include "AliEMCALRecoUtils.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
48 #include "AliGenPythiaEventHeader.h"
49 #include "AliInputEventHandler.h"
50 #include "AliMCEvent.h"
51 #include "AliOADBContainer.h"
52 #include "AliVVertex.h"
53 
54 #include "AliAnalysisTaskChargedParticlesRefMC.h"
55 #include "AliEMCalTriggerWeightHandler.h"
56 
58 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC)
60 
61 namespace EMCalTriggerPtAnalysis {
62 
63 AliAnalysisTaskChargedParticlesRefMC::AliAnalysisTaskChargedParticlesRefMC():
64  AliAnalysisTaskEmcal(),
65  fTrackCuts(nullptr),
66  fTriggerSelection(nullptr),
67  fHistos(nullptr),
68  fWeightHandler(nullptr),
69  fEventTriggers(),
70  fEventWeight(1.),
71  fYshift(0.465),
72  fEtaSign(1),
73  fEtaLabCut(-0.5, 0.5),
74  fEtaCmsCut(-0.13, 0.13),
75  fPhiCut(0., TMath::TwoPi()),
76  fFracPtHard(-1),
77  fNameAcceptanceOADB()
78 {
79  SetCaloTriggerPatchInfoName("EmcalTriggers");
80  SetNeedEmcalGeom(true);
81  SetMakeGeneralHistograms(kTRUE);
82 }
83 
84 AliAnalysisTaskChargedParticlesRefMC::AliAnalysisTaskChargedParticlesRefMC(const char* name):
85  AliAnalysisTaskEmcal(name, true),
86  fTrackCuts(nullptr),
87  fTriggerSelection(nullptr),
88  fHistos(nullptr),
89  fWeightHandler(nullptr),
90  fEventTriggers(),
91  fEventWeight(1.),
92  fYshift(0.465),
93  fEtaSign(1),
94  fEtaLabCut(-0.5, 0.5),
95  fEtaCmsCut(-0.13, 0.13),
96  fPhiCut(0., TMath::TwoPi()),
97  fFracPtHard(-1),
98  fNameAcceptanceOADB()
99 {
100  SetCaloTriggerPatchInfoName("EmcalTriggers");
101  SetNeedEmcalGeom(true);
102  SetMakeGeneralHistograms(kTRUE);
103 }
104 
105 AliAnalysisTaskChargedParticlesRefMC::~AliAnalysisTaskChargedParticlesRefMC() {
106  //if(fTrackCuts) delete fTrackCuts;
107  if(fTriggerSelection) delete fTriggerSelection;
108  if(fHistos) delete fHistos;
109 }
110 
111 void AliAnalysisTaskChargedParticlesRefMC::UserCreateOutputObjects() {
112  AliAnalysisTaskEmcal::UserCreateOutputObjects();
113 
114  if(!fAliAnalysisUtils) fAliAnalysisUtils = new AliAnalysisUtils;
115  fHistos = new THistManager("Ref");
116 
117  if(!fTrackCuts) InitializeTrackCuts("standard",fInputHandler->IsA() == AliAODInputHandler::Class());
118 
119  PtBinning newbinning;
120 
121  fHistos->CreateTH1("hPtHard", "Pt of the hard interaction", 1000, 0., 500);
122  TString triggers[7] = {"True", "MB", "EMC7", "EJ1", "EJ2", "EG1", "EG2"};
123  Double_t ptcuts[5] = {1., 2., 5., 10., 20.};
124  TString species[6] = {"El", "Mu", "Pi", "Ka", "Pr", "Ot"};
125  for(TString *trg = triggers; trg < triggers + sizeof(triggers)/sizeof(TString); trg++){
126  fHistos->CreateTH1(Form("hEventCount%s", trg->Data()), Form("Event Counter for trigger class %s", trg->Data()), 1, 0.5, 1.5);
127  fHistos->CreateTH1(Form("hVertexBefore%s", trg->Data()), Form("Vertex distribution before z-cut for trigger class %s", trg->Data()), 500, -50, 50);
128  fHistos->CreateTH1(Form("hVertexAfter%s", trg->Data()), Form("Vertex distribution after z-cut for trigger class %s", trg->Data()), 100, -10, 10);
129  fHistos->CreateTH1(Form("hPtEtaAll%s", trg->Data()), Form("Charged particle pt distribution all eta trigger %s", trg->Data()), newbinning);
130  fHistos->CreateTH1(Form("hPtEtaCent%s", trg->Data()), Form("Charged particle pt distribution central eta trigger %s", trg->Data()), newbinning);
131  fHistos->CreateTH1(Form("hPtEMCALEtaAll%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta trigger %s", trg->Data()), newbinning);
132  fHistos->CreateTH1(Form("hPtEMCALEtaCent%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta trigger %s", trg->Data()), newbinning);
133  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAll%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution all eta trigger %s", trg->Data()), newbinning);
134  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCent%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution central eta trigger %s", trg->Data()), newbinning);
135  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAll%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution all eta trigger %s", trg->Data()), newbinning);
136  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCent%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution central eta trigger %s", trg->Data()), newbinning);
137  for(TString *piditer = species; piditer < species + sizeof(species)/sizeof(TString); ++piditer){
138  fHistos->CreateTH1(Form("hPtEtaAll%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution all eta trigger %s", piditer->Data(), trg->Data()), newbinning);
139  fHistos->CreateTH1(Form("hPtEtaCent%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution central eta trigger %s", piditer->Data(), trg->Data()), newbinning);
140  fHistos->CreateTH1(Form("hPtEMCALEtaAll%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution all eta trigger %s", piditer->Data(), trg->Data()), newbinning);
141  fHistos->CreateTH1(Form("hPtEMCALEtaCent%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution central eta trigger %s", piditer->Data(), trg->Data()), newbinning);
142  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAll%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution all eta trigger %s", piditer->Data(), trg->Data()), newbinning);
143  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCent%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution central eta trigger %s", piditer->Data(), trg->Data()), newbinning);
144  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAll%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution all eta trigger %s", piditer->Data(), trg->Data()), newbinning);
145  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCent%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution central eta trigger %s", piditer->Data(), trg->Data()), newbinning);
146  }
147  for(int ipt = 0; ipt < 5; ipt++){
148  fHistos->CreateTH1(
149  Form("hEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
150  Form("Eta (lab) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
151  100,
152  -1.,
153  1.
154  );
155  fHistos->CreateTH1(
156  Form("hEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
157  Form("Eta (lab) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
158  100,
159  -1.,
160  1.
161  );
162  fHistos->CreateTH1(
163  Form("hEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
164  Form("Eta (cent) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
165  160,
166  -1.3,
167  1.3
168  );
169  fHistos->CreateTH1(
170  Form("hEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
171  Form("Eta (cent) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
172  160,
173  -1.3,
174  1.3
175  );
176  fHistos->CreateTH1(
177  Form("hEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
178  Form("Eta (lab) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
179  100,
180  -1.,
181  1.
182  );
183  fHistos->CreateTH1(
184  Form("hEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
185  Form("Eta (lab) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
186  100,
187  -1.,
188  1.
189  );
190  fHistos->CreateTH1(
191  Form("hEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
192  Form("Eta (cent) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
193  160,
194  -1.3,
195  1.3
196  );
197  fHistos->CreateTH1(
198  Form("hEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
199  Form("Eta (cent) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
200  160,
201  -1.3,
202  1.3
203  );
204  fHistos->CreateTH1(
205  Form("hPhiDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
206  Form("#phi distribution of particles with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
207  300,
208  0.,
209  2*TMath::Pi()
210  );
211  }
212  }
213  //fHistos->GetListOfHistograms()->Add(fTrackCuts);
214  for(auto hist : *(fHistos->GetListOfHistograms())){
215  fOutput->Add(hist);
216  }
217 
218  PostData(1, fOutput);
219 }
220 
221 bool AliAnalysisTaskChargedParticlesRefMC::IsEventSelected(){
222  fEventTriggers.clear();
223  AliDebugStream(1) << GetName() << ": Using custom event selection" << std::endl;
224  if(!MCEvent()) return false;
225  if(!fTriggerPatchInfo) return false;
226  fEventWeight = fWeightHandler ? fWeightHandler->GetEventWeight(fPythiaHeader) : 1.;
227 
228  // Do MC outlier cut
229  if(fIsPythia){
230  if(!CheckMCOutliers()) return false;
231  }
232 
233 
234  // select trigger
235  bool isMinBias;
236  if((isMinBias = fInputHandler->IsEventSelected() & AliVEvent::kINT7)) fEventTriggers.push_back("MB");
237  // In simulations triggered events are a subset of min. bias events
238  if(fTriggerPatchInfo && fTriggerSelection){
239  if(isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEL0, fTriggerPatchInfo))
240  fEventTriggers.push_back("EMC7"); // triggerstring.Contains("EMC7"),
241  if(isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEJ1, fTriggerPatchInfo))
242  fEventTriggers.push_back("EJ1"); // triggerstring.Contains("EJ1"),
243  if(isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEJ2, fTriggerPatchInfo))
244  fEventTriggers.push_back("EJ2"); // triggerstring.Contains("EJ2"),
245  if(isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEG1, fTriggerPatchInfo))
246  fEventTriggers.push_back("EG1"); // triggerstring.Contains("EG1"),
247  if(isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEG2, fTriggerPatchInfo))
248  fEventTriggers.push_back("EG2"); // triggerstring.Contains("EG2");
249  }
250  if(!fEventTriggers.size()){
251  AliDebugStream(1) << GetName() << ": No trigger selected" << std::endl;
252  }
253 
254  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
255  if(vtx->GetNContributors() < 1) return false;
256  if(!fAliAnalysisUtils->IsVertexSelected2013pA(fInputEvent)) return false; // Apply new vertex cut
257  if(fAliAnalysisUtils->IsPileUpEvent(fInputEvent)) return false; // Apply new vertex cut
258  // Fill reference distribution for the primary vertex before any z-cut
259  fHistos->FillTH1("hVertexBeforeTrue", vtx->GetZ(), fEventWeight);
260  for(const auto &trg : fEventTriggers) fHistos->FillTH1(Form("hVertexBefore%s", trg.c_str()), vtx->GetZ(), fEventWeight);
261  // Apply vertex z cut
262  if(vtx->GetZ() < -10. || vtx->GetZ() > 10.) return false;
263 
264 
265  // Fill Event counter and reference vertex distributions for the different trigger classes
266  fHistos->FillTH1("hEventCountTrue", 1, fEventWeight);
267  fHistos->FillTH1("hVertexAfterTrue", vtx->GetZ(), fEventWeight);
268  for(const auto &trg : fEventTriggers){
269  fHistos->FillTH1(Form("hEventCount%s", trg.c_str()), 1, fEventWeight);
270  fHistos->FillTH1(Form("hVertexAfter%s", trg.c_str()), vtx->GetZ(), fEventWeight);
271  }
272 
273  return true;
274 }
275 
276 void AliAnalysisTaskChargedParticlesRefMC::ExecOnce(){
277  AliAnalysisTaskEmcal::ExecOnce();
278 
279  if(!fLocalInitialized) {
280  AliErrorStream() << GetName() << ": Failed initializing AliAnalysisTaskEmcal" << std::endl;
281  return;
282  }
283 
284  // Load acceptance OADB
285  if(fNameAcceptanceOADB.Length() && fTriggerSelection){
286  AliDebugStream(1) << GetName() << ": Loading acceptance map from OADB file " << fNameAcceptanceOADB << std::endl;
287  AliOADBContainer acceptanceCont("AliEmcalTriggerAcceptance");
288  acceptanceCont.InitFromFile(fNameAcceptanceOADB.Data(), "AliEmcalTriggerAcceptance");
289  TObjArray *acceptanceMaps = dynamic_cast<TObjArray *>(acceptanceCont.GetObject(fInputEvent->GetRunNumber()));
290  TH2 *map(nullptr);
291  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("EG1")))){
292  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger EG1" << std::endl;
293  map->SetDirectory(nullptr);
294  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgEG1, map);
295  }
296  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("EG2")))){
297  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger EG2" << std::endl;
298  map->SetDirectory(nullptr);
299  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgEG2, map);
300  }
301  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("DG1")))){
302  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger DG1" << std::endl;
303  map->SetDirectory(nullptr);
304  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgDG1, map);
305  }
306  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("DG2")))){
307  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger DG2" << std::endl;
308  map->SetDirectory(nullptr);
309  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgDG1, map);
310  }
311  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("EJ1")))){
312  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger EJ1" << std::endl;
313  map->SetDirectory(nullptr);
314  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgEJ1, map);
315  }
316  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("EJ2")))){
317  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger EJ2" << std::endl;
318  map->SetDirectory(nullptr);
319  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgEJ2, map);
320  }
321  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("DJ1")))){
322  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger DJ1" << std::endl;
323  map->SetDirectory(nullptr);
324  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgDJ1, map);
325  }
326  if((map = dynamic_cast<TH2 *>(acceptanceMaps->FindObject("DJ2")))){
327  AliDebugStream(1) << GetName() << ": Found acceptance map for trigger DJ2" << std::endl;
328  map->SetDirectory(nullptr);
329  fTriggerSelection->SetAcceptanceMap(AliEmcalTriggerOfflineSelection::kTrgDJ1, map);
330  }
331  }
332 }
333 
334 bool AliAnalysisTaskChargedParticlesRefMC::Run() {
335  // MonteCarlo Loop
336  // Histograms
337  // - Full eta (-0.8, 0.8), new binning
338  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
339  // - Central eta_{cms} (-0.3, 0.3), new binning,
340  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
341  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
342  AliVParticle *truepart = NULL;
343  Bool_t isEMCAL(kFALSE);
344  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
345  truepart = fMCEvent->GetTrack(ipart);
346 
347  // Select only particles within ALICE acceptance
348  if(!fEtaLabCut.IsInRange(truepart->Eta())) continue;
349  if(!fPhiCut.IsInRange(truepart->Phi())) continue;
350  if(TMath::Abs(truepart->Pt()) < 0.1) continue;
351  if(!truepart->Charge()) continue;
352 
353  if(!IsPhysicalPrimary(truepart, fMCEvent)) continue;
354  isEMCAL = (truepart->Phi() > 1.5 && truepart->Phi() < 3.1) ? kTRUE : kFALSE;
355 
356  // Calculate eta in cms frame according
357  // EPJC74 (2014) 3054:
358  // eta_cms = - eta_lab - |yshift|
359  Double_t etacent = -1. * truepart->Eta() - TMath::Abs(fYshift);
360  etacent *= fEtaSign;
361 
362  Bool_t etacentcut = fEtaCmsCut.IsInRange(etacent);
363 
364  // Get PID
365  TString pid = "";
366  switch(TMath::Abs(truepart->PdgCode())){
367  case kPiPlus: pid = "Pi"; break;
368  case kMuonMinus: pid = "Mu"; break;
369  case kElectron: pid = "El"; break;
370  case kKPlus: pid = "Ka"; break;
371  case kProton: pid = "Pr"; break;
372  default: pid = "Ot"; break;
373  };
374 
375  // Particle selected (do not filter TRD sectors for MC truth)
376  FillTrackHistos("True", fEventWeight, truepart->Pt(), truepart->Eta() * fEtaSign, etacent, truepart->Phi(), etacentcut, isEMCAL, kFALSE, pid);
377  }
378 
379  // Loop over tracks, fill select particles
380  // Histograms
381  // - Full eta (-0.8, 0.8), new binning
382  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
383  // - Central eta (-0.8, -0.2), new binning,
384  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
385  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
386  AliVTrack *checktrack(NULL);
387  AliVParticle *assocMC(NULL);
388  double ptparticle(-1.), etaparticle(-100.), etaEMCAL(0.), phiEMCAL(0.);
389  Bool_t hasTRD = kFALSE;
390  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
391  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
392  if(!checktrack) continue;
393  // Find associated particle
394  assocMC = fMCEvent->GetTrack(TMath::Abs(checktrack->GetLabel()));
395  if(!assocMC) continue; // Fake track
396  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
397 
398  // Select only particles within ALICE acceptance
399  if(!fEtaLabCut.IsInRange(checktrack->Eta())) continue;
400  if(!fPhiCut.IsInRange(checktrack->Phi())) continue;
401  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
402  if(checktrack->IsA() == AliESDtrack::Class()){
403  AliESDtrack copytrack(*(static_cast<AliESDtrack *>(checktrack)));
404  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
405  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
406  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
407  } else {
408  AliAODTrack copytrack(*(static_cast<AliAODTrack *>(checktrack)));
409  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
410  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
411  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
412  }
413  Int_t supermoduleID = -1;
414  isEMCAL = fGeom->SuperModuleNumberFromEtaPhi(etaEMCAL, phiEMCAL, supermoduleID);
415  // Exclude supermodules 10 and 11 as they did not participate in the trigger
416  isEMCAL = isEMCAL && supermoduleID < 10;
417  hasTRD = isEMCAL && supermoduleID >= 4; // supermodules 4 - 10 have TRD in front in the 2012-2013 ALICE setup
418 
419  if(!fTrackCuts->IsTrackAccepted(checktrack)) continue;
420 
421  ptparticle = TMath::Abs(assocMC->Pt());
422  etaparticle = assocMC->Eta();
423 
424  // Calculate eta in cms frame according
425  // EPJC74 (2014) 3054:
426  // eta_cms = - eta_lab - |yshift|
427  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
428  etacent *= fEtaSign;
429 
430  Bool_t etacentcut = fEtaCmsCut.IsInRange(etacent);
431 
432  // Get PID
433  TString assocpid = "";
434  switch(TMath::Abs(assocMC->PdgCode())){
435  case kPiPlus: assocpid = "Pi"; break;
436  case kMuonMinus: assocpid = "Mu"; break;
437  case kElectron: assocpid = "El"; break;
438  case kKPlus: assocpid = "Ka"; break;
439  case kProton: assocpid = "Pr"; break;
440  default: assocpid = "Ot"; break;
441  };
442  for(const auto &trg : fEventTriggers)
443  FillTrackHistos(trg.c_str(), fEventWeight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
444  }
445  return true;
446 }
447 
448 void AliAnalysisTaskChargedParticlesRefMC::FillTrackHistos(
449  const char *eventclass,
450  Double_t weight,
451  Double_t pt,
452  Double_t etalab,
453  Double_t etacent,
454  Double_t phi,
455  Bool_t etacut,
456  Bool_t inEmcal,
457  Bool_t hasTRD,
458  const char *pid
459  )
460 {
461  fHistos->FillTH1(Form("hPtEtaAll%s", eventclass), TMath::Abs(pt), weight);
462  fHistos->FillTH1(Form("hPtEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
463  if(inEmcal){
464  fHistos->FillTH1(Form("hPtEMCALEtaAll%s", eventclass), TMath::Abs(pt), weight);
465  fHistos->FillTH1(Form("hPtEMCALEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
466  if(hasTRD){
467  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAll%s", eventclass), TMath::Abs(pt), weight);
468  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
469  } else {
470  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAll%s", eventclass), TMath::Abs(pt), weight);
471  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAll%s%s", pid, eventclass), TMath::Abs(pt), weight);
472  }
473  }
474 
475  int ptmin[5] = {1,2,5,10,20}; // for eta distributions
476  for(int icut = 0; icut < 5; icut++){
477  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
478  fHistos->FillTH1(Form("hPhiDistAllPt%d%s", ptmin[icut], eventclass), phi, weight);
479  fHistos->FillTH1(Form("hEtaLabDistAllPt%d%s", ptmin[icut], eventclass), etalab, weight);
480  fHistos->FillTH1(Form("hEtaCentDistAllPt%d%s", ptmin[icut], eventclass), etacent, weight);
481  if(inEmcal){
482  fHistos->FillTH1(Form("hEtaLabDistAllEMCALPt%d%s", ptmin[icut], eventclass), etalab, weight);
483  fHistos->FillTH1(Form("hEtaCentDistAllEMCALPt%d%s", ptmin[icut], eventclass), etacent, weight);
484  }
485  }
486  }
487 
488  if(etacut){
489  fHistos->FillTH1(Form("hPtEtaCent%s", eventclass), TMath::Abs(pt), weight);
490  fHistos->FillTH1(Form("hPtEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
491  if(inEmcal){
492  fHistos->FillTH1(Form("hPtEMCALEtaCent%s", eventclass), TMath::Abs(pt), weight);
493  fHistos->FillTH1(Form("hPtEMCALEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
494  if(hasTRD){
495  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCent%s", eventclass), TMath::Abs(pt), weight);
496  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
497  } else {
498  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCent%s", eventclass), TMath::Abs(pt), weight);
499  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCent%s%s", pid, eventclass), TMath::Abs(pt), weight);
500  }
501  }
502  for(int icut = 0; icut < 5; icut++){
503  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
504  fHistos->FillTH1(Form("hEtaLabDistCutPt%d%s", ptmin[icut], eventclass), etalab, weight);
505  fHistos->FillTH1(Form("hEtaCentDistCutPt%d%s", ptmin[icut], eventclass), etacent, weight);
506  if(inEmcal){
507  fHistos->FillTH1(Form("hEtaLabDistCutEMCALPt%d%s", ptmin[icut], eventclass), etalab, weight);
508  fHistos->FillTH1(Form("hEtaCentDistCutEMCALPt%d%s", ptmin[icut], eventclass), etacent, weight);
509  }
510  }
511  }
512  }
513 }
514 
515 
516 void AliAnalysisTaskChargedParticlesRefMC::InitializeTrackCuts(TString cutname, bool isAOD){
517  SetTrackSelection(AliEmcalAnalysisFactory::TrackCutsFactory(cutname, isAOD));
518 }
519 
520 
521 
522 TString AliAnalysisTaskChargedParticlesRefMC::GetFiredTriggerClasses(const TClonesArray* triggerpatches) {
523  TString triggerstring = "";
524  Int_t nEJ1 = 0, nEJ2 = 0, nEG1 = 0, nEG2 = 0;
525  double minADC_EJ1 = 260.,
526  minADC_EJ2 = 127.,
527  minADC_EG1 = 140.,
528  minADC_EG2 = 89.;
529  for(TIter patchIter = TIter(triggerpatches).Begin(); patchIter != TIter::End(); ++patchIter){
530  AliEMCALTriggerPatchInfo *patch = dynamic_cast<AliEMCALTriggerPatchInfo *>(*patchIter);
531  if(!patch->IsOfflineSimple()) continue;
532  if(patch->IsJetHighSimple() && patch->GetADCOfflineAmp() > minADC_EJ1) nEJ1++;
533  if(patch->IsJetLowSimple() && patch->GetADCOfflineAmp() > minADC_EJ2) nEJ2++;
534  if(patch->IsGammaHighSimple() && patch->GetADCOfflineAmp() > minADC_EG1) nEG1++;
535  if(patch->IsGammaLowSimple() && patch->GetADCOfflineAmp() > minADC_EG2) nEG2++;
536  }
537  if(nEJ1) triggerstring += "EJ1";
538  if(nEJ2){
539  if(triggerstring.Length()) triggerstring += ",";
540  triggerstring += "EJ2";
541  }
542  if(nEG1){
543  if(triggerstring.Length()) triggerstring += ",";
544  triggerstring += "EG1";
545  }
546  if(nEG2){
547  if(triggerstring.Length()) triggerstring += ",";
548  triggerstring += "EG2";
549  }
550  return triggerstring;
551 }
552 
553 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) {
554  Bool_t physprim = false;
555  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
556  if(aodmc){
557  physprim = aodmc->IsPhysicalPrimary();
558  } else {
559  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
560  }
561  return physprim;
562 }
563 
567 AliAnalysisTaskChargedParticlesRefMC::PtBinning::PtBinning() :
568  TCustomBinning()
569 {
570  this->SetMinimum(0.);
571  this->AddStep(1, 0.05);
572  this->AddStep(2, 0.1);
573  this->AddStep(4, 0.2);
574  this->AddStep(7, 0.5);
575  this->AddStep(16, 1);
576  this->AddStep(36, 2);
577  this->AddStep(40, 4);
578  this->AddStep(50, 5);
579  this->AddStep(100, 10);
580  this->AddStep(200, 20);
581 }
582 
583 } /* namespace EMCalTriggerPtAnalysis */
AliEMCalTriggerWeightHandler.h
AliAnalysisTaskEmcal::fIsPythia
Bool_t fIsPythia
trigger, if it is a PYTHIA production
Definition: AliAnalysisTaskEmcal.h:300
Double_t
double Double_t
Definition: External.C:58
EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory::TrackCutsFactory
static AliEmcalTrackSelection * TrackCutsFactory(TString name, Bool_t isAOD)
AliEmcalTriggerOfflineSelection.h
AliAnalysisTaskEmcal
Base task in the EMCAL framework.
Definition: AliAnalysisTaskEmcal.h:67
AliAnalysisTaskEmcal::fLocalInitialized
Bool_t fLocalInitialized
whether or not the task has been already initialized
Definition: AliAnalysisTaskEmcal.h:273
THistManager.h
AliAnalysisTaskChargedParticlesRefMC.h
AliEmcalAnalysisFactory.h
AliAnalysisTaskEmcal::SetCaloTriggerPatchInfoName
void SetCaloTriggerPatchInfoName(const char *n)
Definition: AliAnalysisTaskEmcal.h:148
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskEmcal::fGeom
AliEMCALGeometry * fGeom
!emcal geometry
Definition: AliAnalysisTaskEmcal.h:322
AliAnalysisTaskEmcal::fPythiaHeader
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
Definition: AliAnalysisTaskEmcal.h:338
AliAnalysisTaskEmcal::fAliAnalysisUtils
AliAnalysisUtils * fAliAnalysisUtils
!vertex selection (optional)
Definition: AliAnalysisTaskEmcal.h:320
TCustomBinning
Helper class creating user defined custom binning.
Definition: TCustomBinning.h:52
ptmin
const Double_t ptmin
Definition: AddTaskCFDStar.C:38
kProton
Definition: FastCentEstimators.C:23
nullptr
#define nullptr
Definition: AliAnalysisTaskEmcalOfflinePatchesRef.cxx:35
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC
Test class for charged particle distributions (MC case)
Definition: AliAnalysisTaskChargedParticlesRefMC.h:43
AliAnalysisTaskEmcal::IsEventSelected
virtual Bool_t IsEventSelected()
Definition: AliAnalysisTaskEmcal.cxx:1015
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:346
TH2
Definition: External.C:220
ClassImp
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Definition: AliAnalysisTaskCRC.cxx:44
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:159
AliAnalysisTaskEmcal::fTriggerPatchInfo
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Definition: AliAnalysisTaskEmcal.h:327
AliAnalysisTaskEmcal::SetNeedEmcalGeom
void SetNeedEmcalGeom(Bool_t n)
Definition: AliAnalysisTaskEmcal.h:166
THistManager
Container class for histograms for the high- charged particle analysis.
Definition: THistManager.h:43
Bool_t
bool Bool_t
Definition: External.C:53
AliEmcalTrackSelection.h