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