AliPhysics  vAN-20150924 (e816f45)
 All Classes Namespaces Files Functions Variables 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"
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 }
80 
86  AliAnalysisTaskSE(name),
87  fTrackCuts(NULL),
88  fAnalysisUtil(NULL),
89  fHistos(NULL),
90  fGeometry(NULL),
91  fPtHard(0),
92  fPtHardBin(0),
93  fNTrials(0),
94  fXsection(0),
95  fYshift(0.465),
96  fEtaSign(1),
97  fFracPtHard(-1)
98 {
99  // Restrict analysis to the EMCAL acceptance
100  fEtaLabCut[0] = -0.6;
101  fEtaLabCut[1] = 0.6;
102  fEtaCmsCut[0] = -0.13;
103  fEtaCmsCut[1] = 0.13;
104  DefineOutput(1, TList::Class());
105 }
106 
111  if(fTrackCuts) delete fTrackCuts;
112  if(fAnalysisUtil) delete fAnalysisUtil;
113  if(fHistos) delete fHistos;
114 }
119  fAnalysisUtil = new AliAnalysisUtils;
120 
121  fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(true, 1);
122  fTrackCuts->SetName("Standard Track cuts");
123  fTrackCuts->SetMinNCrossedRowsTPC(120);
124  fTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
125 
126  TArrayD oldbinning, newbinning;
127  CreateOldPtBinning(oldbinning);
128  CreateNewPtBinning(newbinning);
129 
130  fHistos = new AliEMCalHistoContainer("Ref");
131  fHistos->CreateTH1("hNtrials", "Number of trials", 1, 0.5, 1.5);
132  fHistos->CreateTProfile("hCrossSection", "PYTHIA cross section", 1, 0.5, 1.5);
133  fHistos->CreateTH1("hNtrialsNoSelect", "Number of trials (without event selection)", 1, 0.5, 1.5);
134  fHistos->CreateTH1("hNtrialsEvent", "Number of trials (from header, after event selection)", 1, 0.5, 1.5);
135  fHistos->CreateTProfile("hCrossSectionEvent", "PYTHIA cross section (from header, after event selection)", 1, 0.5, 1.5);
136  fHistos->CreateTH1("hPtHard", "Pt of the hard interaction", 1000, 0., 500);
137  fHistos->CreateTH1("hTriggerJetPtNoCut", "pt of trigger jets wihtout cuts", 1000, 0., 500);
138  fHistos->CreateTH1("hRatioPtJetPtHardNoCut", "Ratio of pt jet / pt hard without cut", 100, 0., 2.);
139  fHistos->CreateTH1("hTriggerJetPtWithCut", "pt of trigger jets after cuts", 1000, 0., 500);
140  fHistos->CreateTH1("hRatioPtJetPtHardWithCut", "Ratio of pt jet / pt hard with cut on this ratio", 100, 0., 2.);
141  TString triggers[15] = {"True", "MB", "EJ1", "EJ2", "EG1", "EG2", "MBexcl", "EJ2excl", "EG2excl", "E1combined", "E1Jonly", "E1Gonly", "E2combined", "E2Jonly", "E2Gonly"};
142  Double_t ptcuts[5] = {1., 2., 5., 10., 20.};
143  for(TString *trg = triggers; trg < triggers + sizeof(triggers)/sizeof(TString); trg++){
144  fHistos->CreateTH1(Form("hEventCount%s", trg->Data()), Form("Event Counter for trigger class %s", trg->Data()), 1, 0.5, 1.5);
145  fHistos->CreateTH1(Form("hVertexBefore%s", trg->Data()), Form("Vertex distribution before z-cut for trigger class %s", trg->Data()), 500, -50, 50);
146  fHistos->CreateTH1(Form("hVertexAfter%s", trg->Data()), Form("Vertex distribution after z-cut for trigger class %s", trg->Data()), 100, -10, 10);
147  fHistos->CreateTH1(Form("hPtEtaAllOldBinning%s", trg->Data()), Form("Charged particle pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
148  fHistos->CreateTH1(Form("hPtEtaCentOldBinning%s", trg->Data()), Form("Charged particle pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
149  fHistos->CreateTH1(Form("hPtEtaAllNewBinning%s", trg->Data()), Form("Charged particle pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
150  fHistos->CreateTH1(Form("hPtEtaCentNewBinning%s", trg->Data()), Form("Charged particle pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
151  fHistos->CreateTH1(Form("hPtEMCALEtaAllOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
152  fHistos->CreateTH1(Form("hPtEMCALEtaCentOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
153  fHistos->CreateTH1(Form("hPtEMCALEtaAllNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
154  fHistos->CreateTH1(Form("hPtEMCALEtaCentNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
155  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);
156  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);
157  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);
158  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);
159  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);
160  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);
161  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);
162  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);
163  for(int ipt = 0; ipt < 5; ipt++){
165  Form("hEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
166  Form("Eta (lab) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
167  100,
168  -1.,
169  1.
170  );
172  Form("hEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
173  Form("Eta (lab) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
174  100,
175  -1.,
176  1.
177  );
179  Form("hEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
180  Form("Eta (cent) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
181  160,
182  -1.3,
183  1.3
184  );
186  Form("hEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
187  Form("Eta (cent) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
188  160,
189  -1.3,
190  1.3
191  );
193  Form("hEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
194  Form("Eta (lab) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
195  100,
196  -1.,
197  1.
198  );
200  Form("hEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
201  Form("Eta (lab) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
202  100,
203  -1.,
204  1.
205  );
207  Form("hEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
208  Form("Eta (cent) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
209  160,
210  -1.3,
211  1.3
212  );
214  Form("hEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
215  Form("Eta (cent) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
216  160,
217  -1.3,
218  1.3
219  );
221  Form("hPhiDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
222  Form("#phi distribution of particles with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
223  300,
224  0.,
225  2*TMath::Pi()
226  );
227  }
228  }
229  PostData(1, fHistos->GetListOfHistograms());
230 }
231 
232 void AliAnalysisTaskChargedParticlesRefMC::UserExec(Option_t*) { // Select event
233  if(!fMCEvent) return;
234  if(!fGeometry){
235  fGeometry = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
236  }
237  AliGenPythiaEventHeader *pyheader = GetPythiaHeader();
238  if(pyheader){
239  FillTriggerJetHistograms(kFALSE, pyheader);
240  fHistos->FillTH1("hNtrialsNoSelect",1,pyheader->Trials());
241  if(fFracPtHard > 0){
242  // Apply outlier cut in case of a positive fraction of pthard
243  if(IsOutlier(pyheader))
244  return;
245  }
246  FillTriggerJetHistograms(kTRUE, pyheader);
247  }
248  TClonesArray *fTriggerPatches = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject("EmcalTriggers"));
249  if(!fTriggerPatches) return;
250 
251  TString triggerstring = GetFiredTriggerClasses(fTriggerPatches);
252  Bool_t isMinBias = fInputHandler->IsEventSelected() & AliVEvent::kINT7,
253  isEJ1 = isMinBias && triggerstring.Contains("EJ1"),
254  isEJ2 = isMinBias && triggerstring.Contains("EJ2"),
255  isEG1 = isMinBias && triggerstring.Contains("EG1"),
256  isEG2 = isMinBias && triggerstring.Contains("EG2"); // In simulations triggered events are a subset of min. bias events
257  if(!(isMinBias || isEG1 || isEG2 || isEJ1 || isEJ2)) return;
258  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
259  //if(!fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return; // reject pileup event
260  if(vtx->GetNContributors() < 1) return;
261  // Fill reference distribution for the primary vertex before any z-cut
262  fHistos->FillTH1("hVertexBeforeTrue", vtx->GetZ());
263  if(isMinBias) fHistos->FillTH1("hVertexBeforeMB", vtx->GetZ());
264  if(isEJ1) fHistos->FillTH1("hVertexBeforeEJ1", vtx->GetZ());
265  if(isEJ2) fHistos->FillTH1("hVertexBeforeEJ2", vtx->GetZ());
266  if(isEG1) fHistos->FillTH1("hVertexBeforeEG1", vtx->GetZ());
267  if(isEG2) fHistos->FillTH1("hVertexBeforeEG2", vtx->GetZ());
268  if(!fAnalysisUtil->IsVertexSelected2013pA(fInputEvent)) return; // Apply new vertex cut
269  if(fAnalysisUtil->IsPileUpEvent(fInputEvent)) return; // Apply new vertex cut
270  // Apply vertex z cut
271  if(vtx->GetZ() < -10. || vtx->GetZ() > 10.) return;
272 
273  // Fill Event counter and reference vertex distributions for the different trigger classes
274  fHistos->FillTH1("hEventCountTrue", 1);
275  fHistos->FillTH1("hVertexAfterTrue", vtx->GetZ());
276  if(isMinBias){
277  fHistos->FillTH1("hEventCountMB", 1);
278  fHistos->FillTH1("hVertexAfterMB", vtx->GetZ());
279  // check for exclusive classes
280  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
281  fHistos->FillTH1("hEventCountMBexcl", 1);
282  fHistos->FillTH1("hVertexAfterMBexcl", vtx->GetZ());
283  }
284  }
285  if(isEJ1){
286  fHistos->FillTH1("hEventCountEJ1", 1);
287  fHistos->FillTH1("hVertexAfterEJ1", vtx->GetZ());
288  if(isEG1 || isEG2){
289  fHistos->FillTH1("hEventCountE1combined", 1);
290  fHistos->FillTH1("hVertexAfterE1combined", vtx->GetZ());
291  } else {
292  fHistos->FillTH1("hEventCountE1Jonly", 1);
293  fHistos->FillTH1("hVertexAfterE1Jonly", vtx->GetZ());
294  }
295  }
296  if(isEJ2){
297  fHistos->FillTH1("hEventCountEJ2", 1);
298  fHistos->FillTH1("hVertexAfterEJ2", vtx->GetZ());
299  // Check for exclusive classes
300  if(!isEJ1){
301  fHistos->FillTH1("hEventCountEJ2excl", 1);
302  fHistos->FillTH1("hVertexAfterEJ2excl", vtx->GetZ());
303  }
304  if(isEG1 || isEG2){
305  fHistos->FillTH1("hEventCountE2combined", 1);
306  fHistos->FillTH1("hVertexAfterE2combined", vtx->GetZ());
307  } else {
308  fHistos->FillTH1("hEventCountE2Jonly", 1);
309  fHistos->FillTH1("hVertexAfterE2Jonly", vtx->GetZ());
310  }
311  }
312  if(isEG1){
313  fHistos->FillTH1("hEventCountEG1", 1);
314  fHistos->FillTH1("hVertexAfterEG1", vtx->GetZ());
315  if(!(isEJ1 || isEJ2)){
316  fHistos->FillTH1("hEventCountE1Gonly", 1);
317  fHistos->FillTH1("hVertexAfterE1Gonly", vtx->GetZ());
318  }
319  }
320  if(isEG2){
321  fHistos->FillTH1("hEventCountEG2", 1);
322  fHistos->FillTH1("hVertexAfterEG2", vtx->GetZ());
323  // check for exclusive trigger classes
324  if(!isEG1){
325  fHistos->FillTH1("hEventCountEG2excl", 1);
326  fHistos->FillTH1("hVertexAfterEG2excl", vtx->GetZ());
327  }
328  if(!(isEJ1 || isEJ2)){
329  fHistos->FillTH1("hEventCountE2Gonly", 1);
330  fHistos->FillTH1("hVertexAfterE2Gonly", vtx->GetZ());
331  }
332  }
333 
334  // Fill PYTHIA histograms from event header
335  if(pyheader){
336  fHistos->FillTH1("hNtrialsEvent", 1., pyheader->Trials());
337  fHistos->FillProfile("hCrossSectionEvent", 1., pyheader->GetXsection());
338  fHistos->FillTH1("hPtHard", pyheader->GetPtHard());
339  }
340 
341  // MonteCarlo Loop
342  // Histograms
343  // - Full eta (-0.8, 0.8), new binning
344  // - Full eta_{lab} (-0.8, 0.8), old binning
345  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
346  // - Central eta_{cms} (-0.3, 0.3), new binning,
347  // - Central eta_{cms} (-0.3, 0.3), old binning,
348  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
349  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
350  AliVParticle *truepart = NULL;
351  Bool_t isEMCAL(kFALSE);
352  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
353  truepart = fMCEvent->GetTrack(ipart);
354 
355  // Select only particles within ALICE acceptance
356  if((truepart->Eta() < fEtaLabCut[0]) || (truepart->Eta() > fEtaLabCut[1])) continue;
357  if(TMath::Abs(truepart->Pt()) < 0.1) continue;
358  if(!truepart->Charge()) continue;
359 
360  if(!IsPhysicalPrimary(truepart, fMCEvent)) continue;
361  isEMCAL = (truepart->Phi() > 1.5 && truepart->Phi() < 3.1) ? kTRUE : kFALSE;
362 
363  // Calculate eta in cms frame according
364  // EPJC74 (2014) 3054:
365  // eta_cms = - eta_lab - |yshift|
366  Double_t etacent = -1. * truepart->Eta() - TMath::Abs(fYshift);
367  etacent *= fEtaSign;
368 
369  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
370 
371  // Particle selected (do not filter TRD sectors for MC truth)
372  FillTrackHistos("True", truepart->Pt(), truepart->Eta() * fEtaSign, etacent, truepart->Phi(), etacentcut, isEMCAL, kFALSE);
373  }
374 
375  // Loop over tracks, fill select particles
376  // Histograms
377  // - Full eta (-0.8, 0.8), new binning
378  // - Full eta (-0.8, 0.8), old binning
379  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
380  // - Central eta (-0.8, -0.2), new binning,
381  // - Central eta (-0.8, -0.2), old binning,
382  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
383  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
384  AliVTrack *checktrack(NULL);
385  AliVParticle *assocMC(NULL);
386  double ptparticle(-1.), etaparticle(-100.);
387  Bool_t hasTRD = kFALSE;
388  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
389  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
390  if(!checktrack) continue;
391  // Find associated particle
392  assocMC = fMCEvent->GetTrack(TMath::Abs(checktrack->GetLabel()));
393  if(!assocMC) continue; // Fake track
394  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
395 
396  // Select only particles within ALICE acceptance
397  if((checktrack->Eta() < fEtaLabCut[0]) || (checktrack->Eta() > fEtaLabCut[1])) continue;
398  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
399  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(checktrack);
400  Int_t supermoduleID = -1;
401  isEMCAL = fGeometry->SuperModuleNumberFromEtaPhi(checktrack->GetTrackEtaOnEMCal(), checktrack->GetTrackPhiOnEMCal(), supermoduleID);
402  // Exclude supermodules 10 and 11 as they did not participate in the trigger
403  isEMCAL = isEMCAL && supermoduleID < 10;
404  hasTRD = isEMCAL && supermoduleID >= 4; // supermodules 4 - 10 have TRD in front in the 2012-2013 ALICE setup
405 
406  // Distinguish track selection for ESD and AOD tracks
407  AliESDtrack *esdtrack(NULL);
408  AliAODTrack *aodtrack(NULL);
409  if((esdtrack = dynamic_cast<AliESDtrack *>(checktrack))){
410  if(!TrackSelectionESD(esdtrack)) continue;
411  } else if((aodtrack = dynamic_cast<AliAODTrack *>(checktrack))){
412  if(!TrackSelectionAOD(aodtrack)) continue;
413  } else {
414  continue;
415  }
416 
417  ptparticle = TMath::Abs(assocMC->Pt());
418  etaparticle = assocMC->Eta();
419 
420  // Calculate eta in cms frame according
421  // EPJC74 (2014) 3054:
422  // eta_cms = - eta_lab - |yshift|
423  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
424  etacent *= fEtaSign;
425 
426  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
427 
428  // Go through the trigger classes and fill histograms
429  if(isMinBias){
430  FillTrackHistos("MB", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
431  // Check for exclusive classes
432  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
433  FillTrackHistos("MBexcl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
434  }
435  }
436  if(isEJ1){
437  FillTrackHistos("EJ1", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
438  if(isEG1 || isEG2) {
439  FillTrackHistos("E1combined", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
440  } else {
441  FillTrackHistos("E1Jonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
442  }
443  }
444  if(isEJ2){
445  FillTrackHistos("EJ2", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
446  // check for exclusive classes
447  if(!isEJ1){
448  FillTrackHistos("EJ2excl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
449  }
450  if(isEG1 || isEG2) {
451  FillTrackHistos("E2combined", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
452  } else {
453  FillTrackHistos("E2Jonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
454  }
455  }
456  if(isEG1){
457  FillTrackHistos("EG1", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
458  if(!(isEJ1 || isEJ2)){
459  FillTrackHistos("E1Gonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
460  }
461  }
462  if(isEG2){
463  FillTrackHistos("EG2", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
464  // check for exclusive classes
465  if(!isEG1){
466  FillTrackHistos("EG2excl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
467  }
468  if(!(isEJ1 || isEJ2)){
469  FillTrackHistos("E2Gonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
470  }
471  }
472  }
473  PostData(1, fHistos->GetListOfHistograms());
474 
475 }
476 
488  const char *eventclass,
489  Double_t pt,
490  Double_t etalab,
491  Double_t etacent,
492  Double_t phi,
493  Bool_t etacut,
494  Bool_t inEmcal,
495  Bool_t hasTRD
496  )
497 {
498  fHistos->FillTH1(Form("hPtEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
499  fHistos->FillTH1(Form("hPtEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
500  if(inEmcal){
501  fHistos->FillTH1(Form("hPtEMCALEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
502  fHistos->FillTH1(Form("hPtEMCALEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
503  if(hasTRD){
504  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
505  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
506  } else {
507  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
508  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
509  }
510  }
511 
512  int ptmin[5] = {1,2,5,10,20}; // for eta distributions
513  for(int icut = 0; icut < 5; icut++){
514  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
515  fHistos->FillTH1(Form("hPhiDistAllPt%d%s", ptmin[icut], eventclass), phi);
516  fHistos->FillTH1(Form("hEtaLabDistAllPt%d%s", ptmin[icut], eventclass), etalab);
517  fHistos->FillTH1(Form("hEtaCentDistAllPt%d%s", ptmin[icut], eventclass), etacent);
518  if(inEmcal){
519  fHistos->FillTH1(Form("hEtaLabDistAllEMCALPt%d%s", ptmin[icut], eventclass), etalab);
520  fHistos->FillTH1(Form("hEtaCentDistAllEMCALPt%d%s", ptmin[icut], eventclass), etacent);
521  }
522  }
523  }
524 
525  if(etacut){
526  fHistos->FillTH1(Form("hPtEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
527  fHistos->FillTH1(Form("hPtEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
528  if(inEmcal){
529  fHistos->FillTH1(Form("hPtEMCALEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
530  fHistos->FillTH1(Form("hPtEMCALEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
531  if(hasTRD){
532  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
533  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
534  } else {
535  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
536  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
537  }
538  }
539  for(int icut = 0; icut < 5; icut++){
540  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
541  fHistos->FillTH1(Form("hEtaLabDistCutPt%d%s", ptmin[icut], eventclass), etalab);
542  fHistos->FillTH1(Form("hEtaCentDistCutPt%d%s", ptmin[icut], eventclass), etacent);
543  if(inEmcal){
544  fHistos->FillTH1(Form("hEtaLabDistCutEMCALPt%d%s", ptmin[icut], eventclass), etalab);
545  fHistos->FillTH1(Form("hEtaCentDistCutEMCALPt%d%s", ptmin[icut], eventclass), etacent);
546  }
547  }
548  }
549  }
550 }
551 
557 void AliAnalysisTaskChargedParticlesRefMC::FillTriggerJetHistograms(Bool_t aftercut, AliGenPythiaEventHeader *const header){
558  TString ending = aftercut ? "WithCut" : "NoCut";
559  Float_t pbuf[4];
560  TLorentzVector jetvec;
561  TString hnameSpec = "hTriggerJetPt" + ending,
562  hnameptratio = "hRatioPtJetPtHard" + ending;
563  for(int ijet = 0; ijet < header->NTriggerJets(); ijet++){
564  memset(pbuf, 0, sizeof(Float_t) * 4);
565  header->TriggerJet(ijet, pbuf);
566  jetvec.SetPxPyPzE(pbuf[0], pbuf[1], pbuf[2], pbuf[3]);
567  fHistos->FillTH1(hnameSpec.Data(), TMath::Abs(jetvec.Pt()));
568  fHistos->FillTH1(hnameptratio.Data(), TMath::Abs(jetvec.Pt()/header->GetPtHard()));
569  }
570 }
571 
577 {
578  // Called when file changes.
579 
580  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
581  if (!tree) {
582  AliError(Form("%s - UserNotify: No current tree!",GetName()));
583  return kFALSE;
584  }
585 
586  Float_t xsection = 0;
587  Float_t trials = 0;
588  Int_t pthardbin = 0;
589 
590  TFile *curfile = tree->GetCurrentFile();
591  if (!curfile) {
592  AliError(Form("%s - UserNotify: No current file!",GetName()));
593  return kFALSE;
594  }
595 
596  TChain *chain = dynamic_cast<TChain*>(tree);
597  if (chain) tree = chain->GetTree();
598 
599  Int_t nevents = tree->GetEntriesFast();
600 
601  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
602 
603  fHistos->FillTH1("hNtrials", 1., trials);
604  fHistos->FillProfile("hCrossSection", 1., xsection);
605  //fHistos->FillTH1(pthardbin, nevents);
606 
607  return kTRUE;
608 }
609 
622 Bool_t AliAnalysisTaskChargedParticlesRefMC::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard) const {
623 
624  TString file(currFile);
625  fXsec = 0;
626  fTrials = 1;
627 
628  if (file.Contains(".zip#")) {
629  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
630  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
631  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
632  file.Replace(pos+1,pos2-pos1,"");
633  } else {
634  // not an archive take the basename....
635  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
636  }
637  AliDebug(1,Form("File name: %s",file.Data()));
638 
639  // Get the pt hard bin
640  TString strPthard(file);
641 
642  strPthard.Remove(strPthard.Last('/'));
643  strPthard.Remove(strPthard.Last('/'));
644  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
645  strPthard.Remove(0,strPthard.Last('/')+1);
646  if (strPthard.IsDec())
647  pthard = strPthard.Atoi();
648  else
649  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
650 
651  // 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
652  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
653 
654  if (!fxsec) {
655  // next trial fetch the histgram file
656  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
657  if (!fxsec) {
658  // not a severe condition but inciate that we have no information
659  return kFALSE;
660  } else {
661  // find the tlist we want to be independtent of the name so use the Tkey
662  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
663  if (!key) {
664  fxsec->Close();
665  return kFALSE;
666  }
667  TList *list = dynamic_cast<TList*>(key->ReadObj());
668  if (!list) {
669  fxsec->Close();
670  return kFALSE;
671  }
672  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
673  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
674  fxsec->Close();
675  }
676  } else { // no tree pyxsec.root
677  TTree *xtree = (TTree*)fxsec->Get("Xsection");
678  if (!xtree) {
679  fxsec->Close();
680  return kFALSE;
681  }
682  UInt_t ntrials = 0;
683  Double_t xsection = 0;
684  xtree->SetBranchAddress("xsection",&xsection);
685  xtree->SetBranchAddress("ntrials",&ntrials);
686  xtree->GetEntry(0);
687  fTrials = ntrials;
688  fXsec = xsection;
689  fxsec->Close();
690  }
691  return kTRUE;
692 }
693 
699  AliGenPythiaEventHeader *pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
700  if (!pythiaHeader) {
701  // Check if AOD
702  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
703  if (aodMCH) {
704  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
705  pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
706  if (pythiaHeader) break;
707  }
708  }
709  }
710  return pythiaHeader;
711 }
712 
718  std::vector<double> mybinning;
719  std::map<double,double> definitions;
720  definitions.insert(std::pair<double,double>(2.5, 0.1));
721  definitions.insert(std::pair<double,double>(7., 0.25));
722  definitions.insert(std::pair<double,double>(15., 0.5));
723  definitions.insert(std::pair<double,double>(25., 1.));
724  definitions.insert(std::pair<double,double>(40., 2.5));
725  definitions.insert(std::pair<double,double>(50., 5.));
726  definitions.insert(std::pair<double,double>(100., 10.));
727  double currentval = 0;
728  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
729  double limit = id->first, binwidth = id->second;
730  while(currentval < limit){
731  currentval += binwidth;
732  mybinning.push_back(currentval);
733  }
734  }
735  binning.Set(mybinning.size());
736  int ib = 0;
737  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
738  binning[ib++] = *it;
739 }
740 
746  std::vector<double> mybinning;
747  std::map<double,double> definitions;
748  definitions.insert(std::pair<double, double>(1, 0.05));
749  definitions.insert(std::pair<double, double>(2, 0.1));
750  definitions.insert(std::pair<double, double>(4, 0.2));
751  definitions.insert(std::pair<double, double>(7, 0.5));
752  definitions.insert(std::pair<double, double>(16, 1));
753  definitions.insert(std::pair<double, double>(36, 2));
754  definitions.insert(std::pair<double, double>(40, 4));
755  definitions.insert(std::pair<double, double>(50, 5));
756  definitions.insert(std::pair<double, double>(100, 10));
757  definitions.insert(std::pair<double, double>(200, 20));
758  double currentval = 0.;
759  mybinning.push_back(currentval);
760  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
761  double limit = id->first, binwidth = id->second;
762  while(currentval < limit){
763  currentval += binwidth;
764  mybinning.push_back(currentval);
765  }
766  }
767  binning.Set(mybinning.size());
768  int ib = 0;
769  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
770  binning[ib++] = *it;
771 }
772 
779  return fTrackCuts->AcceptTrack(track);
780 }
781 
788  if(!track->TestFilterBit(AliAODTrack::kTrkGlobal)) return false;
789  if(track->GetTPCNCrossedRows() < 120) return false;
790  return true;
791 }
792 
798 TString AliAnalysisTaskChargedParticlesRefMC::GetFiredTriggerClasses(const TClonesArray* triggerpatches) {
799  TString triggerstring = "";
800  Int_t nEJ1 = 0, nEJ2 = 0, nEG1 = 0, nEG2 = 0;
801  double minADC_EJ1 = 260.,
802  minADC_EJ2 = 127.,
803  minADC_EG1 = 140.,
804  minADC_EG2 = 89.;
805  for(TIter patchIter = TIter(triggerpatches).Begin(); patchIter != TIter::End(); ++patchIter){
806  AliEmcalTriggerPatchInfo *patch = dynamic_cast<AliEmcalTriggerPatchInfo *>(*patchIter);
807  if(!patch->IsOfflineSimple()) continue;
808  if(patch->IsJetHighSimple() && patch->GetADCOfflineAmp() > minADC_EJ1) nEJ1++;
809  if(patch->IsJetLowSimple() && patch->GetADCOfflineAmp() > minADC_EJ2) nEJ2++;
810  if(patch->IsGammaHighSimple() && patch->GetADCOfflineAmp() > minADC_EG1) nEG1++;
811  if(patch->IsGammaLowSimple() && patch->GetADCOfflineAmp() > minADC_EG2) nEG2++;
812  }
813  if(nEJ1) triggerstring += "EJ1";
814  if(nEJ2){
815  if(triggerstring.Length()) triggerstring += ",";
816  triggerstring += "EJ2";
817  }
818  if(nEG1){
819  if(triggerstring.Length()) triggerstring += ",";
820  triggerstring += "EG1";
821  }
822  if(nEG2){
823  if(triggerstring.Length()) triggerstring += ",";
824  triggerstring += "EG2";
825  }
826  return triggerstring;
827 }
828 
837 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) {
838  Bool_t physprim = false;
839  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
840  if(aodmc){
841  physprim = aodmc->IsPhysicalPrimary();
842  } else {
843  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
844  }
845  return physprim;
846 }
847 
853 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsOutlier(AliGenPythiaEventHeader * const header) const {
854  Bool_t hasOutlier = kFALSE;
855  Float_t pbuf[4];
856  TLorentzVector jetvec;
857  for(int ijet = 0; ijet < header->NTriggerJets(); ijet++){
858  memset(pbuf, 0, sizeof(Float_t) * 4);
859  header->TriggerJet(ijet, pbuf);
860  jetvec.SetPxPyPzE(pbuf[0], pbuf[1], pbuf[2], pbuf[3]);
861  if(TMath::Abs(jetvec.Pt()) >= this->fFracPtHard * header->GetPtHard()){
862  hasOutlier = true;
863  break;
864  }
865  }
866  return hasOutlier;
867 }
868 
869 } /* namespace EMCalTriggerPtAnalysis */
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
void FillTriggerJetHistograms(Bool_t aftercut, AliGenPythiaEventHeader *const header)
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.
Main data structure storing all relevant information of EMCAL/DCAL trigger patches.
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.
Class to make array of trigger patch objects in AOD/ESD events.
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.)