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