AliPhysics  c69cd47 (c69cd47)
 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 <vector>
16 #include <map>
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 "AliVVertex.h"
52 
53 #include "AliAnalysisTaskChargedParticlesRefMC.h"
54 #include "AliEMCalTriggerWeightHandler.h"
55 
57 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC)
59 
60 namespace EMCalTriggerPtAnalysis {
61 
65 AliAnalysisTaskChargedParticlesRefMC::AliAnalysisTaskChargedParticlesRefMC():
66  AliAnalysisTaskSE(),
67  fTrackCuts(NULL),
68  fAnalysisUtil(NULL),
69  fTriggerSelection(NULL),
70  fHistos(NULL),
71  fGeometry(NULL),
72  fWeightHandler(NULL),
73  fUsePythiaHard(kFALSE),
74  fPtHard(0),
75  fPtHardBin(0),
76  fNTrials(0),
77  fXsection(0),
78  fYshift(0.465),
79  fEtaSign(1),
80  fFracPtHard(-1)
81 {
82  // Restrict analysis to the EMCAL acceptance
83  fEtaLabCut[0] = -0.6;
84  fEtaLabCut[1] = 0.6;
85  fEtaCmsCut[0] = -0.13;
86  fEtaCmsCut[1] = 0.13;
87 }
88 
93 AliAnalysisTaskChargedParticlesRefMC::AliAnalysisTaskChargedParticlesRefMC(const char* name):
94  AliAnalysisTaskSE(name),
95  fTrackCuts(NULL),
96  fAnalysisUtil(NULL),
97  fTriggerSelection(NULL),
98  fHistos(NULL),
99  fGeometry(NULL),
100  fWeightHandler(NULL),
101  fUsePythiaHard(kFALSE),
102  fPtHard(0),
103  fPtHardBin(0),
104  fNTrials(0),
105  fXsection(0),
106  fYshift(0.465),
107  fEtaSign(1),
108  fFracPtHard(-1)
109 {
110  // Restrict analysis to the EMCAL acceptance
111  fEtaLabCut[0] = -0.6;
112  fEtaLabCut[1] = 0.6;
113  fEtaCmsCut[0] = -0.13;
114  fEtaCmsCut[1] = 0.13;
115  DefineOutput(1, TList::Class());
116 }
117 
121 AliAnalysisTaskChargedParticlesRefMC::~AliAnalysisTaskChargedParticlesRefMC() {
122  //if(fTrackCuts) delete fTrackCuts;
123  if(fAnalysisUtil) delete fAnalysisUtil;
124  if(fTriggerSelection) delete fTriggerSelection;
125  if(fHistos) delete fHistos;
126 }
130 void AliAnalysisTaskChargedParticlesRefMC::UserCreateOutputObjects() {
131  if(!fAnalysisUtil) fAnalysisUtil = new AliAnalysisUtils;
132  fHistos = new THistManager("Ref");
133 
134  if(!fTrackCuts) InitializeTrackCuts("standard",fInputHandler->IsA() == AliAODInputHandler::Class());
135 
136  TArrayD oldbinning, newbinning;
137  CreateOldPtBinning(oldbinning);
138  CreateNewPtBinning(newbinning);
139 
140  fHistos->CreateTH1("hNtrials", "Number of trials", 1, 0.5, 1.5);
141  fHistos->CreateTProfile("hCrossSection", "PYTHIA cross section", 1, 0.5, 1.5);
142  fHistos->CreateTH1("hNtrialsNoSelect", "Number of trials (without event selection)", 1, 0.5, 1.5);
143  fHistos->CreateTH1("hNtrialsEvent", "Number of trials (from header, after event selection)", 1, 0.5, 1.5);
144  fHistos->CreateTProfile("hCrossSectionEvent", "PYTHIA cross section (from header, after event selection)", 1, 0.5, 1.5);
145  fHistos->CreateTH1("hPtHard", "Pt of the hard interaction", 1000, 0., 500);
146  fHistos->CreateTH1("hTriggerJetPtNoCut", "pt of trigger jets wihtout cuts", 1000, 0., 500);
147  fHistos->CreateTH1("hRatioPtJetPtHardNoCut", "Ratio of pt jet / pt hard without cut", 1000, 0., 20.);
148  fHistos->CreateTH1("hTriggerJetPtWithCut", "pt of trigger jets after cuts", 1000, 0., 500);
149  fHistos->CreateTH1("hRatioPtJetPtHardWithCut", "Ratio of pt jet / pt hard with cut on this ratio", 1000, 0., 20.);
150  TString triggers[16] = {"True", "MB", "EMC7", "EJ1", "EJ2", "EG1", "EG2", "MBexcl", "EJ2excl", "EG2excl", "E1combined", "E1Jonly", "E1Gonly", "E2combined", "E2Jonly", "E2Gonly"};
151  Double_t ptcuts[5] = {1., 2., 5., 10., 20.};
152  TString species[6] = {"El", "Mu", "Pi", "Ka", "Pr", "Ot"};
153  for(TString *trg = triggers; trg < triggers + sizeof(triggers)/sizeof(TString); trg++){
154  fHistos->CreateTH1(Form("hEventCount%s", trg->Data()), Form("Event Counter for trigger class %s", trg->Data()), 1, 0.5, 1.5);
155  fHistos->CreateTH1(Form("hVertexBefore%s", trg->Data()), Form("Vertex distribution before z-cut for trigger class %s", trg->Data()), 500, -50, 50);
156  fHistos->CreateTH1(Form("hVertexAfter%s", trg->Data()), Form("Vertex distribution after z-cut for trigger class %s", trg->Data()), 100, -10, 10);
157  fHistos->CreateTH1(Form("hPtEtaAllOldBinning%s", trg->Data()), Form("Charged particle pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
158  fHistos->CreateTH1(Form("hPtEtaCentOldBinning%s", trg->Data()), Form("Charged particle pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
159  fHistos->CreateTH1(Form("hPtEtaAllNewBinning%s", trg->Data()), Form("Charged particle pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
160  fHistos->CreateTH1(Form("hPtEtaCentNewBinning%s", trg->Data()), Form("Charged particle pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
161  fHistos->CreateTH1(Form("hPtEMCALEtaAllOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
162  fHistos->CreateTH1(Form("hPtEMCALEtaCentOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
163  fHistos->CreateTH1(Form("hPtEMCALEtaAllNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
164  fHistos->CreateTH1(Form("hPtEMCALEtaCentNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
165  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAllOldBinning%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
166  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCentOldBinning%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
167  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAllNewBinning%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
168  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCentNewBinning%s", trg->Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
169  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAllOldBinning%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
170  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCentOldBinning%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
171  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAllNewBinning%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
172  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCentNewBinning%s", trg->Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
173  for(TString *piditer = species; piditer < species + sizeof(species)/sizeof(TString); ++piditer){
174  fHistos->CreateTH1(Form("hPtEtaAllOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution all eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
175  fHistos->CreateTH1(Form("hPtEtaCentOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution central eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
176  fHistos->CreateTH1(Form("hPtEtaAllNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution all eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
177  fHistos->CreateTH1(Form("hPtEtaCentNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s pt distribution central eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
178  fHistos->CreateTH1(Form("hPtEMCALEtaAllOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution all eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
179  fHistos->CreateTH1(Form("hPtEMCALEtaCentOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution central eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
180  fHistos->CreateTH1(Form("hPtEMCALEtaAllNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution all eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
181  fHistos->CreateTH1(Form("hPtEMCALEtaCentNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL pt distribution central eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
182  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAllOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution all eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
183  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCentOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution central eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
184  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAllNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution all eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
185  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCentNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (no TRD in front) pt distribution central eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
186  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAllOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution all eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
187  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCentOldBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution central eta old binning trigger %s", piditer->Data(), trg->Data()), oldbinning);
188  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAllNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution all eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
189  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCentNewBinning%s%s", piditer->Data(), trg->Data()), Form("Charged %s in EMCAL (with TRD in front) pt distribution central eta new binning trigger %s", piditer->Data(), trg->Data()), newbinning);
190  }
191  for(int ipt = 0; ipt < 5; ipt++){
192  fHistos->CreateTH1(
193  Form("hEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
194  Form("Eta (lab) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
195  100,
196  -1.,
197  1.
198  );
199  fHistos->CreateTH1(
200  Form("hEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
201  Form("Eta (lab) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
202  100,
203  -1.,
204  1.
205  );
206  fHistos->CreateTH1(
207  Form("hEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
208  Form("Eta (cent) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
209  160,
210  -1.3,
211  1.3
212  );
213  fHistos->CreateTH1(
214  Form("hEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
215  Form("Eta (cent) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
216  160,
217  -1.3,
218  1.3
219  );
220  fHistos->CreateTH1(
221  Form("hEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
222  Form("Eta (lab) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
223  100,
224  -1.,
225  1.
226  );
227  fHistos->CreateTH1(
228  Form("hEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
229  Form("Eta (lab) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
230  100,
231  -1.,
232  1.
233  );
234  fHistos->CreateTH1(
235  Form("hEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
236  Form("Eta (cent) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
237  160,
238  -1.3,
239  1.3
240  );
241  fHistos->CreateTH1(
242  Form("hEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
243  Form("Eta (cent) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
244  160,
245  -1.3,
246  1.3
247  );
248  fHistos->CreateTH1(
249  Form("hPhiDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
250  Form("#phi distribution of particles with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
251  300,
252  0.,
253  2*TMath::Pi()
254  );
255  }
256  }
257  //fHistos->GetListOfHistograms()->Add(fTrackCuts);
258  PostData(1, fHistos->GetListOfHistograms());
259 }
260 
261 void AliAnalysisTaskChargedParticlesRefMC::UserExec(Option_t*) { // Select event
262  if(!fMCEvent) return;
263  if(!fGeometry){
264  fGeometry = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
265  }
266  AliGenPythiaEventHeader *pyheader = GetPythiaHeader();
267  if(pyheader){
268  FillTriggerJetHistograms(kFALSE, pyheader);
269  fHistos->FillTH1("hNtrialsNoSelect",1,pyheader->Trials());
270  if(fFracPtHard > 0){
271  // Apply outlier cut in case of a positive fraction of pthard
272  if(IsOutlier(pyheader))
273  return;
274  }
275  FillTriggerJetHistograms(kTRUE, pyheader);
276  }
277 
278  Double_t weight = fWeightHandler ? fWeightHandler->GetEventWeight(pyheader) : 1.;
279 
280  //TString triggerstring = GetFiredTriggerClasses(fTriggerPatches);
281  Bool_t isMinBias = fInputHandler->IsEventSelected() & AliVEvent::kINT7,
282  isEMC7 = kFALSE, isEJ1 = kFALSE, isEJ2 = kFALSE, isEG1 = kFALSE, isEG2 = kFALSE;
283  TClonesArray *fTriggerPatches = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject("EmcalTriggers"));
284  if(fTriggerPatches && fTriggerSelection){
285  isEMC7 = isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEL0, fTriggerPatches); // triggerstring.Contains("EMC7"),
286  isEJ1 = isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEJ1, fTriggerPatches); // triggerstring.Contains("EJ1"),
287  isEJ2 = isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEJ2, fTriggerPatches); // triggerstring.Contains("EJ2"),
288  isEG1 = isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEG1, fTriggerPatches); // triggerstring.Contains("EG1"),
289  isEG2 = isMinBias && fTriggerSelection->IsOfflineSelected(AliEmcalTriggerOfflineSelection::kTrgEG2, fTriggerPatches); // triggerstring.Contains("EG2"); // In simulations triggered events are a subset of min. bias events
290 
291  }
292  if(!(isMinBias || isEG1 || isEG2 || isEJ1 || isEJ2)) return;
293  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
294  //if(!fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return; // reject pileup event
295  if(vtx->GetNContributors() < 1) return;
296  // Fill reference distribution for the primary vertex before any z-cut
297  fHistos->FillTH1("hVertexBeforeTrue", vtx->GetZ(), weight);
298  if(isMinBias) fHistos->FillTH1("hVertexBeforeMB", vtx->GetZ(), weight);
299  if(isEMC7) fHistos->FillTH1("hVertexBeforeEMC7", vtx->GetZ(), weight);
300  if(isEJ1) fHistos->FillTH1("hVertexBeforeEJ1", vtx->GetZ(), weight);
301  if(isEJ2) fHistos->FillTH1("hVertexBeforeEJ2", vtx->GetZ(), weight);
302  if(isEG1) fHistos->FillTH1("hVertexBeforeEG1", vtx->GetZ(), weight);
303  if(isEG2) fHistos->FillTH1("hVertexBeforeEG2", vtx->GetZ(), weight);
304  if(!fAnalysisUtil->IsVertexSelected2013pA(fInputEvent)) return; // Apply new vertex cut
305  if(fAnalysisUtil->IsPileUpEvent(fInputEvent)) return; // Apply new vertex cut
306  // Apply vertex z cut
307  if(vtx->GetZ() < -10. || vtx->GetZ() > 10.) return;
308 
309  // Fill Event counter and reference vertex distributions for the different trigger classes
310  fHistos->FillTH1("hEventCountTrue", 1, weight);
311  fHistos->FillTH1("hVertexAfterTrue", vtx->GetZ(), weight);
312  if(isMinBias){
313  fHistos->FillTH1("hEventCountMB", 1, weight);
314  fHistos->FillTH1("hVertexAfterMB", vtx->GetZ(), weight);
315  // check for exclusive classes
316  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
317  fHistos->FillTH1("hEventCountMBexcl", 1, weight);
318  fHistos->FillTH1("hVertexAfterMBexcl", vtx->GetZ(), weight);
319  }
320  }
321  if(isEMC7){
322  fHistos->FillTH1("hEventCountEMC7", 1, weight);
323  fHistos->FillTH1("hVertexAfterEMC7", vtx->GetZ(), weight);
324  }
325  if(isEJ1){
326  fHistos->FillTH1("hEventCountEJ1", 1, weight);
327  fHistos->FillTH1("hVertexAfterEJ1", vtx->GetZ(), weight);
328  if(isEG1 || isEG2){
329  fHistos->FillTH1("hEventCountE1combined", 1, weight);
330  fHistos->FillTH1("hVertexAfterE1combined", vtx->GetZ(), weight);
331  } else {
332  fHistos->FillTH1("hEventCountE1Jonly", 1, weight);
333  fHistos->FillTH1("hVertexAfterE1Jonly", vtx->GetZ(), weight);
334  }
335  }
336  if(isEJ2){
337  fHistos->FillTH1("hEventCountEJ2", 1, weight);
338  fHistos->FillTH1("hVertexAfterEJ2", vtx->GetZ(), weight);
339  // Check for exclusive classes
340  if(!isEJ1){
341  fHistos->FillTH1("hEventCountEJ2excl", 1, weight);
342  fHistos->FillTH1("hVertexAfterEJ2excl", vtx->GetZ(), weight);
343  }
344  if(isEG1 || isEG2){
345  fHistos->FillTH1("hEventCountE2combined", 1, weight);
346  fHistos->FillTH1("hVertexAfterE2combined", vtx->GetZ(), weight);
347  } else {
348  fHistos->FillTH1("hEventCountE2Jonly", 1, weight);
349  fHistos->FillTH1("hVertexAfterE2Jonly", vtx->GetZ(), weight);
350  }
351  }
352  if(isEG1){
353  fHistos->FillTH1("hEventCountEG1", 1, weight);
354  fHistos->FillTH1("hVertexAfterEG1", vtx->GetZ(), weight);
355  if(!(isEJ1 || isEJ2)){
356  fHistos->FillTH1("hEventCountE1Gonly", 1, weight);
357  fHistos->FillTH1("hVertexAfterE1Gonly", vtx->GetZ(), weight);
358  }
359  }
360  if(isEG2){
361  fHistos->FillTH1("hEventCountEG2", 1, weight);
362  fHistos->FillTH1("hVertexAfterEG2", vtx->GetZ(), weight);
363  // check for exclusive trigger classes
364  if(!isEG1){
365  fHistos->FillTH1("hEventCountEG2excl", 1, weight);
366  fHistos->FillTH1("hVertexAfterEG2excl", vtx->GetZ(), weight);
367  }
368  if(!(isEJ1 || isEJ2)){
369  fHistos->FillTH1("hEventCountE2Gonly", 1, weight);
370  fHistos->FillTH1("hVertexAfterE2Gonly", vtx->GetZ(), weight);
371  }
372  }
373 
374  // Fill PYTHIA histograms from event header
375  if(pyheader){
376  fHistos->FillTH1("hNtrialsEvent", 1., pyheader->Trials());
377  fHistos->FillProfile("hCrossSectionEvent", 1., pyheader->GetXsection());
378  fHistos->FillTH1("hPtHard", pyheader->GetPtHard());
379  }
380 
381  // MonteCarlo Loop
382  // Histograms
383  // - Full eta (-0.8, 0.8), new binning
384  // - Full eta_{lab} (-0.8, 0.8), old binning
385  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
386  // - Central eta_{cms} (-0.3, 0.3), new binning,
387  // - Central eta_{cms} (-0.3, 0.3), old binning,
388  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
389  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
390  AliVParticle *truepart = NULL;
391  Bool_t isEMCAL(kFALSE);
392  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
393  truepart = fMCEvent->GetTrack(ipart);
394 
395  // Select only particles within ALICE acceptance
396  if((truepart->Eta() < fEtaLabCut[0]) || (truepart->Eta() > fEtaLabCut[1])) continue;
397  if(TMath::Abs(truepart->Pt()) < 0.1) continue;
398  if(!truepart->Charge()) continue;
399 
400  if(!IsPhysicalPrimary(truepart, fMCEvent)) continue;
401  isEMCAL = (truepart->Phi() > 1.5 && truepart->Phi() < 3.1) ? kTRUE : kFALSE;
402 
403  // Calculate eta in cms frame according
404  // EPJC74 (2014) 3054:
405  // eta_cms = - eta_lab - |yshift|
406  Double_t etacent = -1. * truepart->Eta() - TMath::Abs(fYshift);
407  etacent *= fEtaSign;
408 
409  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
410 
411  // Get PID
412  TString pid = "";
413  switch(TMath::Abs(truepart->PdgCode())){
414  case kPiPlus: pid = "Pi"; break;
415  case kMuonMinus: pid = "Mu"; break;
416  case kElectron: pid = "El"; break;
417  case kKPlus: pid = "Ka"; break;
418  case kProton: pid = "Pr"; break;
419  default: pid = "Ot"; break;
420  };
421 
422  // Particle selected (do not filter TRD sectors for MC truth)
423  FillTrackHistos("True", weight, truepart->Pt(), truepart->Eta() * fEtaSign, etacent, truepart->Phi(), etacentcut, isEMCAL, kFALSE, pid);
424  }
425 
426  // Loop over tracks, fill select particles
427  // Histograms
428  // - Full eta (-0.8, 0.8), new binning
429  // - Full eta (-0.8, 0.8), old binning
430  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
431  // - Central eta (-0.8, -0.2), new binning,
432  // - Central eta (-0.8, -0.2), old binning,
433  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
434  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
435  AliVTrack *checktrack(NULL);
436  AliVParticle *assocMC(NULL);
437  double ptparticle(-1.), etaparticle(-100.), etaEMCAL(0.), phiEMCAL(0.);
438  Bool_t hasTRD = kFALSE;
439  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
440  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
441  if(!checktrack) continue;
442  // Find associated particle
443  assocMC = fMCEvent->GetTrack(TMath::Abs(checktrack->GetLabel()));
444  if(!assocMC) continue; // Fake track
445  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
446 
447  // Select only particles within ALICE acceptance
448  if((checktrack->Eta() < fEtaLabCut[0]) || (checktrack->Eta() > fEtaLabCut[1])) continue;
449  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
450  if(checktrack->IsA() == AliESDtrack::Class()){
451  AliESDtrack copytrack(*(static_cast<AliESDtrack *>(checktrack)));
452  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
453  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
454  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
455  } else {
456  AliAODTrack copytrack(*(static_cast<AliAODTrack *>(checktrack)));
457  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
458  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
459  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
460  }
461  Int_t supermoduleID = -1;
462  isEMCAL = fGeometry->SuperModuleNumberFromEtaPhi(etaEMCAL, phiEMCAL, supermoduleID);
463  // Exclude supermodules 10 and 11 as they did not participate in the trigger
464  isEMCAL = isEMCAL && supermoduleID < 10;
465  hasTRD = isEMCAL && supermoduleID >= 4; // supermodules 4 - 10 have TRD in front in the 2012-2013 ALICE setup
466 
467  if(!fTrackCuts->IsTrackAccepted(checktrack)) continue;
468 
469  ptparticle = TMath::Abs(assocMC->Pt());
470  etaparticle = assocMC->Eta();
471 
472  // Calculate eta in cms frame according
473  // EPJC74 (2014) 3054:
474  // eta_cms = - eta_lab - |yshift|
475  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
476  etacent *= fEtaSign;
477 
478  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
479 
480  // Get PID
481  TString assocpid = "";
482  switch(TMath::Abs(assocMC->PdgCode())){
483  case kPiPlus: assocpid = "Pi"; break;
484  case kMuonMinus: assocpid = "Mu"; break;
485  case kElectron: assocpid = "El"; break;
486  case kKPlus: assocpid = "Ka"; break;
487  case kProton: assocpid = "Pr"; break;
488  default: assocpid = "Ot"; break;
489  };
490 
491  // Go through the trigger classes and fill histograms
492  if(isMinBias){
493  FillTrackHistos("MB", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
494  // Check for exclusive classes
495  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
496  FillTrackHistos("MBexcl", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
497  }
498  }
499  if(isEMC7){
500  FillTrackHistos("EMC7", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
501  }
502  if(isEJ1){
503  FillTrackHistos("EJ1", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
504  if(isEG1 || isEG2) {
505  FillTrackHistos("E1combined", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
506  } else {
507  FillTrackHistos("E1Jonly", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
508  }
509  }
510  if(isEJ2){
511  FillTrackHistos("EJ2", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
512  // check for exclusive classes
513  if(!isEJ1){
514  FillTrackHistos("EJ2excl", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
515  }
516  if(isEG1 || isEG2) {
517  FillTrackHistos("E2combined", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
518  } else {
519  FillTrackHistos("E2Jonly", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
520  }
521  }
522  if(isEG1){
523  FillTrackHistos("EG1", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
524  if(!(isEJ1 || isEJ2)){
525  FillTrackHistos("E1Gonly", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
526  }
527  }
528  if(isEG2){
529  FillTrackHistos("EG2", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
530  // check for exclusive classes
531  if(!isEG1){
532  FillTrackHistos("EG2excl", weight, ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
533  }
534  if(!(isEJ1 || isEJ2)){
535  FillTrackHistos("E2Gonly", weight, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD, assocpid);
536  }
537  }
538  }
539  PostData(1, fHistos->GetListOfHistograms());
540 
541 }
542 
554 void AliAnalysisTaskChargedParticlesRefMC::FillTrackHistos(
555  const char *eventclass,
556  Double_t weight,
557  Double_t pt,
558  Double_t etalab,
559  Double_t etacent,
560  Double_t phi,
561  Bool_t etacut,
562  Bool_t inEmcal,
563  Bool_t hasTRD,
564  const char *pid
565  )
566 {
567  fHistos->FillTH1(Form("hPtEtaAllNewBinning%s", eventclass), TMath::Abs(pt), weight);
568  fHistos->FillTH1(Form("hPtEtaAllOldBinning%s", eventclass), TMath::Abs(pt), weight);
569  fHistos->FillTH1(Form("hPtEtaAllNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
570  fHistos->FillTH1(Form("hPtEtaAllOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
571  if(inEmcal){
572  fHistos->FillTH1(Form("hPtEMCALEtaAllNewBinning%s", eventclass), TMath::Abs(pt), weight);
573  fHistos->FillTH1(Form("hPtEMCALEtaAllOldBinning%s", eventclass), TMath::Abs(pt), weight);
574  fHistos->FillTH1(Form("hPtEMCALEtaAllNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
575  fHistos->FillTH1(Form("hPtEMCALEtaAllOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
576  if(hasTRD){
577  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllNewBinning%s", eventclass), TMath::Abs(pt), weight);
578  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllOldBinning%s", eventclass), TMath::Abs(pt), weight);
579  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
580  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
581  } else {
582  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllNewBinning%s", eventclass), TMath::Abs(pt), weight);
583  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllOldBinning%s", eventclass), TMath::Abs(pt), weight);
584  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
585  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
586  }
587  }
588 
589  int ptmin[5] = {1,2,5,10,20}; // for eta distributions
590  for(int icut = 0; icut < 5; icut++){
591  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
592  fHistos->FillTH1(Form("hPhiDistAllPt%d%s", ptmin[icut], eventclass), phi, weight);
593  fHistos->FillTH1(Form("hEtaLabDistAllPt%d%s", ptmin[icut], eventclass), etalab, weight);
594  fHistos->FillTH1(Form("hEtaCentDistAllPt%d%s", ptmin[icut], eventclass), etacent, weight);
595  if(inEmcal){
596  fHistos->FillTH1(Form("hEtaLabDistAllEMCALPt%d%s", ptmin[icut], eventclass), etalab, weight);
597  fHistos->FillTH1(Form("hEtaCentDistAllEMCALPt%d%s", ptmin[icut], eventclass), etacent, weight);
598  }
599  }
600  }
601 
602  if(etacut){
603  fHistos->FillTH1(Form("hPtEtaCentNewBinning%s", eventclass), TMath::Abs(pt), weight);
604  fHistos->FillTH1(Form("hPtEtaCentOldBinning%s", eventclass), TMath::Abs(pt), weight);
605  fHistos->FillTH1(Form("hPtEtaCentNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
606  fHistos->FillTH1(Form("hPtEtaCentOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
607  if(inEmcal){
608  fHistos->FillTH1(Form("hPtEMCALEtaCentNewBinning%s", eventclass), TMath::Abs(pt), weight);
609  fHistos->FillTH1(Form("hPtEMCALEtaCentOldBinning%s", eventclass), TMath::Abs(pt), weight);
610  fHistos->FillTH1(Form("hPtEMCALEtaCentNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
611  fHistos->FillTH1(Form("hPtEMCALEtaCentOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
612  if(hasTRD){
613  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentNewBinning%s", eventclass), TMath::Abs(pt), weight);
614  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentOldBinning%s", eventclass), TMath::Abs(pt), weight);
615  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
616  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
617  } else {
618  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentNewBinning%s", eventclass), TMath::Abs(pt), weight);
619  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentOldBinning%s", eventclass), TMath::Abs(pt), weight);
620  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentNewBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
621  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentOldBinning%s%s", pid, eventclass), TMath::Abs(pt), weight);
622  }
623  }
624  for(int icut = 0; icut < 5; icut++){
625  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
626  fHistos->FillTH1(Form("hEtaLabDistCutPt%d%s", ptmin[icut], eventclass), etalab, weight);
627  fHistos->FillTH1(Form("hEtaCentDistCutPt%d%s", ptmin[icut], eventclass), etacent, weight);
628  if(inEmcal){
629  fHistos->FillTH1(Form("hEtaLabDistCutEMCALPt%d%s", ptmin[icut], eventclass), etalab, weight);
630  fHistos->FillTH1(Form("hEtaCentDistCutEMCALPt%d%s", ptmin[icut], eventclass), etacent, weight);
631  }
632  }
633  }
634  }
635 }
636 
642 void AliAnalysisTaskChargedParticlesRefMC::FillTriggerJetHistograms(Bool_t aftercut, AliGenPythiaEventHeader *const header){
643  TString ending = aftercut ? "WithCut" : "NoCut";
644  Float_t pbuf[4];
645  TLorentzVector jetvec;
646  TString hnameSpec = "hTriggerJetPt" + ending,
647  hnameptratio = "hRatioPtJetPtHard" + ending;
648  for(int ijet = 0; ijet < header->NTriggerJets(); ijet++){
649  memset(pbuf, 0, sizeof(Float_t) * 4);
650  header->TriggerJet(ijet, pbuf);
651  jetvec.SetPxPyPzE(pbuf[0], pbuf[1], pbuf[2], pbuf[3]);
652  fHistos->FillTH1(hnameSpec.Data(), TMath::Abs(jetvec.Pt()));
653  fHistos->FillTH1(hnameptratio.Data(), TMath::Abs(jetvec.Pt()/header->GetPtHard()));
654  }
655 }
656 
661 Bool_t AliAnalysisTaskChargedParticlesRefMC::UserNotify()
662 {
663  // Called when file changes.
664  if(!fUsePythiaHard) return kTRUE;
665 
666  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
667  if (!tree) {
668  AliError(Form("%s - UserNotify: No current tree!",GetName()));
669  return kFALSE;
670  }
671 
672  Float_t xsection = 0;
673  Float_t trials = 0;
674  Int_t pthardbin = 0;
675 
676  TFile *curfile = tree->GetCurrentFile();
677  if (!curfile) {
678  AliError(Form("%s - UserNotify: No current file!",GetName()));
679  return kFALSE;
680  }
681 
682  TChain *chain = dynamic_cast<TChain*>(tree);
683  if (chain) tree = chain->GetTree();
684 
685  Int_t nevents = tree->GetEntriesFast();
686 
687  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
688 
689  fHistos->FillTH1("hNtrials", 1., trials);
690  fHistos->FillProfile("hCrossSection", 1., xsection);
691  //fHistos->FillTH1(pthardbin, nevents);
692 
693  return kTRUE;
694 }
695 
708 Bool_t AliAnalysisTaskChargedParticlesRefMC::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard) const {
709 
710  TString file(currFile);
711  fXsec = 0;
712  fTrials = 1;
713 
714  if (file.Contains(".zip#")) {
715  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
716  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
717  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
718  file.Replace(pos+1,pos2-pos1,"");
719  } else {
720  // not an archive take the basename....
721  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
722  }
723  AliDebug(1,Form("File name: %s",file.Data()));
724 
725  // Get the pt hard bin
726  TString strPthard(file);
727 
728  strPthard.Remove(strPthard.Last('/'));
729  strPthard.Remove(strPthard.Last('/'));
730  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
731  strPthard.Remove(0,strPthard.Last('/')+1);
732  if (strPthard.IsDec())
733  pthard = strPthard.Atoi();
734  else
735  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
736 
737  // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
738  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
739 
740  if (!fxsec) {
741  // next trial fetch the histgram file
742  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
743  if (!fxsec) {
744  // not a severe condition but inciate that we have no information
745  return kFALSE;
746  } else {
747  // find the tlist we want to be independtent of the name so use the Tkey
748  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
749  if (!key) {
750  fxsec->Close();
751  return kFALSE;
752  }
753  TList *list = dynamic_cast<TList*>(key->ReadObj());
754  if (!list) {
755  fxsec->Close();
756  return kFALSE;
757  }
758  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
759  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
760  fxsec->Close();
761  }
762  } else { // no tree pyxsec.root
763  TTree *xtree = (TTree*)fxsec->Get("Xsection");
764  if (!xtree) {
765  fxsec->Close();
766  return kFALSE;
767  }
768  UInt_t ntrials = 0;
769  Double_t xsection = 0;
770  xtree->SetBranchAddress("xsection",&xsection);
771  xtree->SetBranchAddress("ntrials",&ntrials);
772  xtree->GetEntry(0);
773  fTrials = ntrials;
774  fXsec = xsection;
775  fxsec->Close();
776  }
777  return kTRUE;
778 }
779 
784 AliGenPythiaEventHeader *AliAnalysisTaskChargedParticlesRefMC::GetPythiaHeader() const {
785  AliGenPythiaEventHeader *pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
786  if (!pythiaHeader) {
787  // Check if AOD
788  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
789  if (aodMCH) {
790  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
791  pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
792  if (pythiaHeader) break;
793  }
794  }
795  }
796  return pythiaHeader;
797 }
798 
804 void AliAnalysisTaskChargedParticlesRefMC::InitializeTrackCuts(TString cutname, bool isAOD){
805  SetTrackSelection(AliEmcalAnalysisFactory::TrackCutsFactory(cutname, isAOD));
806 }
807 
812 void AliAnalysisTaskChargedParticlesRefMC::CreateOldPtBinning(TArrayD& binning) const {
813  std::vector<double> mybinning;
814  std::map<double,double> definitions;
815  definitions.insert(std::pair<double,double>(2.5, 0.1));
816  definitions.insert(std::pair<double,double>(7., 0.25));
817  definitions.insert(std::pair<double,double>(15., 0.5));
818  definitions.insert(std::pair<double,double>(25., 1.));
819  definitions.insert(std::pair<double,double>(40., 2.5));
820  definitions.insert(std::pair<double,double>(50., 5.));
821  definitions.insert(std::pair<double,double>(100., 10.));
822  double currentval = 0;
823  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
824  double limit = id->first, binwidth = id->second;
825  while(currentval < limit){
826  currentval += binwidth;
827  mybinning.push_back(currentval);
828  }
829  }
830  binning.Set(mybinning.size());
831  int ib = 0;
832  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
833  binning[ib++] = *it;
834 }
835 
840 void AliAnalysisTaskChargedParticlesRefMC::CreateNewPtBinning(TArrayD& binning) const {
841  std::vector<double> mybinning;
842  std::map<double,double> definitions;
843  definitions.insert(std::pair<double, double>(1, 0.05));
844  definitions.insert(std::pair<double, double>(2, 0.1));
845  definitions.insert(std::pair<double, double>(4, 0.2));
846  definitions.insert(std::pair<double, double>(7, 0.5));
847  definitions.insert(std::pair<double, double>(16, 1));
848  definitions.insert(std::pair<double, double>(36, 2));
849  definitions.insert(std::pair<double, double>(40, 4));
850  definitions.insert(std::pair<double, double>(50, 5));
851  definitions.insert(std::pair<double, double>(100, 10));
852  definitions.insert(std::pair<double, double>(200, 20));
853  double currentval = 0.;
854  mybinning.push_back(currentval);
855  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
856  double limit = id->first, binwidth = id->second;
857  while(currentval < limit){
858  currentval += binwidth;
859  mybinning.push_back(currentval);
860  }
861  }
862  binning.Set(mybinning.size());
863  int ib = 0;
864  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
865  binning[ib++] = *it;
866 }
867 
873 TString AliAnalysisTaskChargedParticlesRefMC::GetFiredTriggerClasses(const TClonesArray* triggerpatches) {
874  TString triggerstring = "";
875  Int_t nEJ1 = 0, nEJ2 = 0, nEG1 = 0, nEG2 = 0;
876  double minADC_EJ1 = 260.,
877  minADC_EJ2 = 127.,
878  minADC_EG1 = 140.,
879  minADC_EG2 = 89.;
880  for(TIter patchIter = TIter(triggerpatches).Begin(); patchIter != TIter::End(); ++patchIter){
881  AliEMCALTriggerPatchInfo *patch = dynamic_cast<AliEMCALTriggerPatchInfo *>(*patchIter);
882  if(!patch->IsOfflineSimple()) continue;
883  if(patch->IsJetHighSimple() && patch->GetADCOfflineAmp() > minADC_EJ1) nEJ1++;
884  if(patch->IsJetLowSimple() && patch->GetADCOfflineAmp() > minADC_EJ2) nEJ2++;
885  if(patch->IsGammaHighSimple() && patch->GetADCOfflineAmp() > minADC_EG1) nEG1++;
886  if(patch->IsGammaLowSimple() && patch->GetADCOfflineAmp() > minADC_EG2) nEG2++;
887  }
888  if(nEJ1) triggerstring += "EJ1";
889  if(nEJ2){
890  if(triggerstring.Length()) triggerstring += ",";
891  triggerstring += "EJ2";
892  }
893  if(nEG1){
894  if(triggerstring.Length()) triggerstring += ",";
895  triggerstring += "EG1";
896  }
897  if(nEG2){
898  if(triggerstring.Length()) triggerstring += ",";
899  triggerstring += "EG2";
900  }
901  return triggerstring;
902 }
903 
912 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) {
913  Bool_t physprim = false;
914  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
915  if(aodmc){
916  physprim = aodmc->IsPhysicalPrimary();
917  } else {
918  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
919  }
920  return physprim;
921 }
922 
928 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsOutlier(AliGenPythiaEventHeader * const header) const {
929  Bool_t hasOutlier = kFALSE;
930  Float_t pbuf[4];
931  TLorentzVector jetvec;
932  for(int ijet = 0; ijet < header->NTriggerJets(); ijet++){
933  memset(pbuf, 0, sizeof(Float_t) * 4);
934  header->TriggerJet(ijet, pbuf);
935  jetvec.SetPxPyPzE(pbuf[0], pbuf[1], pbuf[2], pbuf[3]);
936  if(TMath::Abs(jetvec.Pt()) >= this->fFracPtHard * header->GetPtHard()){
937  hasOutlier = true;
938  break;
939  }
940  }
941  return hasOutlier;
942 }
943 
944 } /* namespace EMCalTriggerPtAnalysis */
AliEMCalTriggerWeightHandler.h
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::fWeightHandler
const AliEMCalTriggerWeightHandler * fWeightHandler
Weight handler (optional)
Definition: AliAnalysisTaskChargedParticlesRefMC.h:83
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::FillTriggerJetHistograms
void FillTriggerJetHistograms(Bool_t aftercut, AliGenPythiaEventHeader *const header)
Definition: AliAnalysisTaskChargedParticlesRefMC.cxx:642
EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory::TrackCutsFactory
static AliEmcalTrackSelection * TrackCutsFactory(TString name, Bool_t isAOD)
AliEmcalTriggerOfflineSelection.h
AliAnalysisTaskChargedParticlesRefMC.h
gSystem
TSystem * gSystem
list
TList * list
Definition: DrawAnaCaloTrackQA.C:18
AliEmcalAnalysisFactory.h
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::fFracPtHard
Double_t fFracPtHard
Cut on the maximum fraction of pt hard of any trigger jet.
Definition: AliAnalysisTaskChargedParticlesRefMC.h:96
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::fTriggerSelection
AliEmcalTriggerOfflineSelection * fTriggerSelection
Offline trigger selection.
Definition: AliAnalysisTaskChargedParticlesRefMC.h:80
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::IsPhysicalPrimary
Bool_t IsPhysicalPrimary(const AliVParticle *const part, AliMCEvent *const mcevent)
Definition: AliAnalysisTaskChargedParticlesRefMC.cxx:912
ptmin
const Double_t ptmin
Definition: AddTaskCFDStar.C:38
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::PythiaInfoFromFile
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard) const
Definition: AliAnalysisTaskChargedParticlesRefMC.cxx:708
EMCalTriggerPtAnalysis::AliEmcalTriggerOfflineSelection::IsOfflineSelected
Bool_t IsOfflineSelected(EmcalTriggerClass trgcls, const TClonesArray *const triggerpatches) const
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC
Unit test class for charged particle distributions (MC case)
Definition: AliAnalysisTaskChargedParticlesRefMC.h:35
EMCalTriggerPtAnalysis::AliEMCalTriggerWeightHandler::GetEventWeight
double GetEventWeight(const AliMCEvent *const event) const
Definition: AliEMCalTriggerWeightHandler.cxx:66
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::FillTrackHistos
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, Bool_t hasTRD, const char *pid)
Definition: AliAnalysisTaskChargedParticlesRefMC.cxx:554
ClassImp
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Definition: AliAnalysisTaskCRC.cxx:44
file
TFile * file
Definition: DrawAnaCaloTrackQA.C:19
nevents
Int_t nevents[nsamples]
Definition: FitMassSpectra.C:74
EMCalTriggerPtAnalysis::AliAnalysisTaskChargedParticlesRefMC::fEtaSign
Double_t fEtaSign
Sign of the eta distribution (swaps when beam directions swap): p-Pb: +1, Pb-p: -1.
Definition: AliAnalysisTaskChargedParticlesRefMC.h:93
GetRunNumber
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
AliEmcalTrackSelection.h
AliEmcalTrackSelection::IsTrackAccepted
virtual bool IsTrackAccepted(AliVTrack *const trk)=0