AliPhysics  v5-06-40-01 (42bb456)
 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 "AliESDtrackCuts.h"
39 #include "AliESDEvent.h"
40 #include "AliESDtrack.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliInputEventHandler.h"
43 #include "AliMCEvent.h"
44 #include "AliVVertex.h"
45 
46 #include "AliEMCalHistoContainer.h"
48 
52 
53 namespace EMCalTriggerPtAnalysis {
54 
58 AliAnalysisTaskChargedParticlesRefMC::AliAnalysisTaskChargedParticlesRefMC():
59  AliAnalysisTaskSE(),
60  fTrackCuts(NULL),
61  fAnalysisUtil(NULL),
62  fHistos(NULL),
63  fPtHard(0),
64  fPtHardBin(0),
65  fNTrials(0),
66  fXsection(0),
67  fYshift(0.465),
68  fEtaSign(1)
69 {
70  // Restrict analysis to the EMCAL acceptance
71  fEtaLabCut[0] = -0.6;
72  fEtaLabCut[1] = 0.6;
73  fEtaCmsCut[0] = -0.13;
74  fEtaCmsCut[1] = 0.13;
75 }
76 
82  AliAnalysisTaskSE(name),
83  fTrackCuts(NULL),
84  fAnalysisUtil(NULL),
85  fHistos(NULL),
86  fPtHard(0),
87  fPtHardBin(0),
88  fNTrials(0),
89  fXsection(0),
90  fYshift(0.465),
91  fEtaSign(1)
92 {
93  // Restrict analysis to the EMCAL acceptance
94  fEtaLabCut[0] = -0.6;
95  fEtaLabCut[1] = 0.6;
96  fEtaCmsCut[0] = -0.13;
97  fEtaCmsCut[1] = 0.13;
98  DefineOutput(1, TList::Class());
99 }
100 
105  if(fTrackCuts) delete fTrackCuts;
106  if(fAnalysisUtil) delete fAnalysisUtil;
107  if(fHistos) delete fHistos;
108 }
113  fAnalysisUtil = new AliAnalysisUtils;
114 
115  fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(true, 1);
116  fTrackCuts->SetName("Standard Track cuts");
117  fTrackCuts->SetMinNCrossedRowsTPC(120);
118  fTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
119 
120  TArrayD oldbinning, newbinning;
121  CreateOldPtBinning(oldbinning);
122  CreateNewPtBinning(newbinning);
123 
124  fHistos = new AliEMCalHistoContainer("Ref");
125  fHistos->CreateTH1("hNtrials", "Number of trials", 1, 0.5, 1.5);
126  fHistos->CreateTProfile("hCrossSection", "PYTHIA cross section", 1, 0.5, 1.5);
127  fHistos->CreateTH1("hNtrialsNoSelect", "Number of trials (without event selection)", 1, 0.5, 1.5);
128  fHistos->CreateTH1("hNtrialsEvent", "Number of trials (from header, after event selection)", 1, 0.5, 1.5);
129  fHistos->CreateTProfile("hCrossSectionEvent", "PYTHIA cross section (from header, after event selection)", 1, 0.5, 1.5);
130  fHistos->CreateTH1("hPtHard", "Pt of the hard interaction", 1000, 0., 500);
131  TString triggers[15] = {"True", "MB", "EJ1", "EJ2", "EG1", "EG2", "MBexcl", "EJ2excl", "EG2excl", "E1combined", "E1Jonly", "E1Gonly", "E2combined", "E2Jonly", "E2Gonly"};
132  Double_t ptcuts[5] = {1., 2., 5., 10., 20.};
133  for(TString *trg = triggers; trg < triggers + sizeof(triggers)/sizeof(TString); trg++){
134  fHistos->CreateTH1(Form("hEventCount%s", trg->Data()), Form("Event Counter for trigger class %s", trg->Data()), 1, 0.5, 1.5);
135  fHistos->CreateTH1(Form("hVertexBefore%s", trg->Data()), Form("Vertex distribution before z-cut for trigger class %s", trg->Data()), 500, -50, 50);
136  fHistos->CreateTH1(Form("hVertexAfter%s", trg->Data()), Form("Vertex distribution after z-cut for trigger class %s", trg->Data()), 100, -10, 10);
137  fHistos->CreateTH1(Form("hPtEtaAllOldBinning%s", trg->Data()), Form("Charged particle pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
138  fHistos->CreateTH1(Form("hPtEtaCentOldBinning%s", trg->Data()), Form("Charged particle pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
139  fHistos->CreateTH1(Form("hPtEtaAllNewBinning%s", trg->Data()), Form("Charged particle pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
140  fHistos->CreateTH1(Form("hPtEtaCentNewBinning%s", trg->Data()), Form("Charged particle pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
141  fHistos->CreateTH1(Form("hPtEMCALEtaAllOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta old binning trigger %s", trg->Data()), oldbinning);
142  fHistos->CreateTH1(Form("hPtEMCALEtaCentOldBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta old binning trigger %s", trg->Data()), oldbinning);
143  fHistos->CreateTH1(Form("hPtEMCALEtaAllNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution all eta new binning trigger %s", trg->Data()), newbinning);
144  fHistos->CreateTH1(Form("hPtEMCALEtaCentNewBinning%s", trg->Data()), Form("Charged particle in EMCAL pt distribution central eta new binning trigger %s", trg->Data()), newbinning);
145  for(int ipt = 0; ipt < 5; ipt++){
147  Form("hEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
148  Form("Eta (lab) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
149  100,
150  -1.,
151  1.
152  );
154  Form("hEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
155  Form("Eta (lab) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
156  100,
157  -1.,
158  1.
159  );
161  Form("hEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
162  Form("Eta (cent) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
163  160,
164  -1.3,
165  1.3
166  );
168  Form("hEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
169  Form("Eta (cent) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
170  160,
171  -1.3,
172  1.3
173  );
175  Form("hEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
176  Form("Eta (lab) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
177  100,
178  -1.,
179  1.
180  );
182  Form("hEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
183  Form("Eta (lab) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
184  100,
185  -1.,
186  1.
187  );
189  Form("hEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
190  Form("Eta (cent) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
191  160,
192  -1.3,
193  1.3
194  );
196  Form("hEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
197  Form("Eta (cent) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
198  160,
199  -1.3,
200  1.3
201  );
203  Form("hPhiDistAllPt%d%s", static_cast<Int_t>(ptcuts[ipt]), trg->Data()),
204  Form("#phi distribution of particles with Pt above %.1f GeV/c trigger %s", ptcuts[ipt], trg->Data()),
205  300,
206  0.,
207  2*TMath::Pi()
208  );
209  }
210  }
211  PostData(1, fHistos->GetListOfHistograms());
212 }
213 
214 void AliAnalysisTaskChargedParticlesRefMC::UserExec(Option_t*) { // Select event
215  if(!fMCEvent) return;
216  AliGenPythiaEventHeader *pyheader = GetPythiaHeader();
217  if(pyheader){
218  fHistos->FillTH1("hNtrialsNoSelect",1,pyheader->Trials());
219  }
220  TClonesArray *fTriggerPatches = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject("EmcalTriggers"));
221  if(!fTriggerPatches) return;
222 
223  TString triggerstring = GetFiredTriggerClasses(fTriggerPatches);
224  Bool_t isMinBias = fInputHandler->IsEventSelected() & AliVEvent::kINT7,
225  isEJ1 = triggerstring.Contains("EJ1"),
226  isEJ2 = triggerstring.Contains("EJ2"),
227  isEG1 = triggerstring.Contains("EG1"),
228  isEG2 = triggerstring.Contains("EG2");
229  if(!(isMinBias || isEG1 || isEG2 || isEJ1 || isEJ2)) return;
230  const AliVVertex *vtx = fInputEvent->GetPrimaryVertex();
231  //if(!fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.)) return; // reject pileup event
232  if(vtx->GetNContributors() < 1) return;
233  // Fill reference distribution for the primary vertex before any z-cut
234  fHistos->FillTH1("hVertexBeforeTrue", vtx->GetZ());
235  if(isMinBias) fHistos->FillTH1("hVertexBeforeMB", vtx->GetZ());
236  if(isEJ1) fHistos->FillTH1("hVertexBeforeEJ1", vtx->GetZ());
237  if(isEJ2) fHistos->FillTH1("hVertexBeforeEJ2", vtx->GetZ());
238  if(isEG1) fHistos->FillTH1("hVertexBeforeEG1", vtx->GetZ());
239  if(isEG2) fHistos->FillTH1("hVertexBeforeEG2", vtx->GetZ());
240  if(!fAnalysisUtil->IsVertexSelected2013pA(fInputEvent)) return; // Apply new vertex cut
241  if(fAnalysisUtil->IsPileUpEvent(fInputEvent)) return; // Apply new vertex cut
242  // Apply vertex z cut
243  if(vtx->GetZ() < -10. || vtx->GetZ() > 10.) return;
244 
245  // Fill Event counter and reference vertex distributions for the different trigger classes
246  fHistos->FillTH1("hEventCountTrue", 1);
247  fHistos->FillTH1("hVertexAfterTrue", vtx->GetZ());
248  if(isMinBias){
249  fHistos->FillTH1("hEventCountMB", 1);
250  fHistos->FillTH1("hVertexAfterMB", vtx->GetZ());
251  // check for exclusive classes
252  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
253  fHistos->FillTH1("hEventCountMBexcl", 1);
254  fHistos->FillTH1("hVertexAfterMBexcl", vtx->GetZ());
255  }
256  }
257  if(isEJ1){
258  fHistos->FillTH1("hEventCountEJ1", 1);
259  fHistos->FillTH1("hVertexAfterEJ1", vtx->GetZ());
260  if(isEG1 || isEG2){
261  fHistos->FillTH1("hEventCountE1combined", 1);
262  fHistos->FillTH1("hVertexAfterE1combined", vtx->GetZ());
263  } else {
264  fHistos->FillTH1("hEventCountE1Jonly", 1);
265  fHistos->FillTH1("hVertexAfterE1Jonly", vtx->GetZ());
266  }
267  }
268  if(isEJ2){
269  fHistos->FillTH1("hEventCountEJ2", 1);
270  fHistos->FillTH1("hVertexAfterEJ2", vtx->GetZ());
271  // Check for exclusive classes
272  if(!isEJ1){
273  fHistos->FillTH1("hEventCountEJ2excl", 1);
274  fHistos->FillTH1("hVertexAfterEJ2excl", vtx->GetZ());
275  }
276  if(isEG1 || isEG2){
277  fHistos->FillTH1("hEventCountE2combined", 1);
278  fHistos->FillTH1("hVertexAfterE2combined", vtx->GetZ());
279  } else {
280  fHistos->FillTH1("hEventCountE2Jonly", 1);
281  fHistos->FillTH1("hVertexAfterE2Jonly", vtx->GetZ());
282  }
283  }
284  if(isEG1){
285  fHistos->FillTH1("hEventCountEG1", 1);
286  fHistos->FillTH1("hVertexAfterEG1", vtx->GetZ());
287  if(!(isEJ1 || isEJ2)){
288  fHistos->FillTH1("hEventCountE1Gonly", 1);
289  fHistos->FillTH1("hVertexAfterE1Gonly", vtx->GetZ());
290  }
291  }
292  if(isEG2){
293  fHistos->FillTH1("hEventCountEG2", 1);
294  fHistos->FillTH1("hVertexAfterEG2", vtx->GetZ());
295  // check for exclusive trigger classes
296  if(!isEG1){
297  fHistos->FillTH1("hEventCountEG2excl", 1);
298  fHistos->FillTH1("hVertexAfterEG2excl", vtx->GetZ());
299  }
300  if(!(isEJ1 || isEJ2)){
301  fHistos->FillTH1("hEventCountE2Gonly", 1);
302  fHistos->FillTH1("hVertexAfterE2Gonly", vtx->GetZ());
303  }
304  }
305 
306  // Fill PYTHIA histograms from event header
307  if(pyheader){
308  fHistos->FillTH1("hNtrialsEvent", 1., pyheader->Trials());
309  fHistos->FillProfile("hCrossSectionEvent", 1., pyheader->GetXsection());
310  fHistos->FillTH1("hPtHard", pyheader->GetPtHard());
311  }
312 
313  // MonteCarlo Loop
314  // Histograms
315  // - Full eta (-0.8, 0.8), new binning
316  // - Full eta_{lab} (-0.8, 0.8), old binning
317  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
318  // - Central eta_{cms} (-0.3, 0.3), new binning,
319  // - Central eta_{cms} (-0.3, 0.3), old binning,
320  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
321  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
322  AliVParticle *truepart = NULL;
323  Bool_t isEMCAL(kFALSE);
324  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
325  truepart = fMCEvent->GetTrack(ipart);
326 
327  // Select only particles within ALICE acceptance
328  if((truepart->Eta() < fEtaLabCut[0]) || (truepart->Eta() > fEtaLabCut[1])) continue;
329  if(TMath::Abs(truepart->Pt()) < 0.1) continue;
330  if(!truepart->Charge()) continue;
331 
332  if(!IsPhysicalPrimary(truepart, fMCEvent)) continue;
333 
334  isEMCAL = (truepart->Phi() > 1.5 && truepart->Phi() < 3.1) ? kTRUE : kFALSE;
335 
336  // Calculate eta in cms frame according
337  // EPJC74 (2014) 3054:
338  // eta_cms = - eta_lab - |yshift|
339  Double_t etacent = -1. * truepart->Eta() - TMath::Abs(fYshift);
340  etacent *= fEtaSign;
341 
342  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
343 
344  // Particle selected
345  FillTrackHistos("True", truepart->Pt(), truepart->Eta() * fEtaSign, etacent, truepart->Phi(), etacentcut, isEMCAL);
346  }
347 
348  // Loop over tracks, fill select particles
349  // Histograms
350  // - Full eta (-0.8, 0.8), new binning
351  // - Full eta (-0.8, 0.8), old binning
352  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
353  // - Central eta (-0.8, -0.2), new binning,
354  // - Central eta (-0.8, -0.2), old binning,
355  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
356  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
357  AliVTrack *checktrack(NULL);
358  AliVParticle *assocMC(NULL);
359  double ptparticle(-1.), etaparticle(-100.);
360  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
361  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
362  if(!checktrack) continue;
363  // Find associated particle
364  assocMC = fMCEvent->GetTrack(TMath::Abs(checktrack->GetLabel()));
365  if(!assocMC) continue; // Fake track
366  if(!IsPhysicalPrimary(assocMC, fMCEvent)) continue;
367 
368  // Select only particles within ALICE acceptance
369  if((checktrack->Eta() < fEtaLabCut[0]) || (checktrack->Eta() > fEtaLabCut[1])) continue;
370  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
371  isEMCAL = (checktrack->Phi() > 1.5 && checktrack->Phi() < 3.1) ? kTRUE : kFALSE;
372 
373  // Distinguish track selection for ESD and AOD tracks
374  AliESDtrack *esdtrack(NULL);
375  AliAODTrack *aodtrack(NULL);
376  if((esdtrack = dynamic_cast<AliESDtrack *>(checktrack))){
377  if(!TrackSelectionESD(esdtrack)) continue;
378  } else if((aodtrack = dynamic_cast<AliAODTrack *>(checktrack))){
379  if(!TrackSelectionAOD(aodtrack)) continue;
380  } else {
381  continue;
382  }
383 
384  ptparticle = TMath::Abs(assocMC->Pt());
385  etaparticle = assocMC->Eta();
386 
387  // Calculate eta in cms frame according
388  // EPJC74 (2014) 3054:
389  // eta_cms = - eta_lab - |yshift|
390  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
391  etacent *= fEtaSign;
392 
393  Bool_t etacentcut = etacent > fEtaCmsCut[0] && etacent < fEtaCmsCut[1];
394 
395  // Go through the trigger classes and fill histograms
396  if(isMinBias){
397  FillTrackHistos("MB", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
398  // Check for exclusive classes
399  if(!(isEG1 || isEG2 || isEJ1 || isEJ2)){
400  FillTrackHistos("MBexcl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
401  }
402  }
403  if(isEJ1){
404  FillTrackHistos("EJ1", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
405  if(isEG1 || isEG2) {
406  FillTrackHistos("E1combined", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
407  } else {
408  FillTrackHistos("E1Jonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
409  }
410  }
411  if(isEJ2){
412  FillTrackHistos("EJ2", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
413  // check for exclusive classes
414  if(!isEJ1){
415  FillTrackHistos("EJ2excl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
416  }
417  if(isEG1 || isEG2) {
418  FillTrackHistos("E2combined", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
419  } else {
420  FillTrackHistos("E2Jonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
421  }
422  }
423  if(isEG1){
424  FillTrackHistos("EG1", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
425  if(!(isEJ1 || isEJ2)){
426  FillTrackHistos("E1Gonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
427  }
428  }
429  if(isEG2){
430  FillTrackHistos("EG2", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
431  // check for exclusive classes
432  if(!isEG1){
433  FillTrackHistos("EG2excl", ptparticle, checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
434  }
435  if(!(isEJ1 || isEJ2)){
436  FillTrackHistos("E2Gonly", checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL);
437  }
438  }
439  }
440  PostData(1, fHistos->GetListOfHistograms());
441 
442 }
443 
455  const char *eventclass,
456  Double_t pt,
457  Double_t etalab,
458  Double_t etacent,
459  Double_t phi,
460  Bool_t etacut,
461  Bool_t inEmcal
462  )
463 {
464  fHistos->FillTH1(Form("hPtEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
465  fHistos->FillTH1(Form("hPtEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
466  if(inEmcal){
467  fHistos->FillTH1(Form("hPtEMCALEtaAllNewBinning%s", eventclass), TMath::Abs(pt));
468  fHistos->FillTH1(Form("hPtEMCALEtaAllOldBinning%s", eventclass), TMath::Abs(pt));
469  }
470 
471  int ptmin[5] = {1,2,5,10,20}; // for eta distributions
472  for(int icut = 0; icut < 5; icut++){
473  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
474  fHistos->FillTH1(Form("hPhiDistAllPt%d%s", ptmin[icut], eventclass), phi);
475  fHistos->FillTH1(Form("hEtaLabDistAllPt%d%s", ptmin[icut], eventclass), etalab);
476  fHistos->FillTH1(Form("hEtaCentDistAllPt%d%s", ptmin[icut], eventclass), etacent);
477  if(inEmcal){
478  fHistos->FillTH1(Form("hEtaLabDistAllEMCALPt%d%s", ptmin[icut], eventclass), etalab);
479  fHistos->FillTH1(Form("hEtaCentDistAllEMCALPt%d%s", ptmin[icut], eventclass), etacent);
480  }
481  }
482  }
483 
484  if(etacut){
485  fHistos->FillTH1(Form("hPtEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
486  fHistos->FillTH1(Form("hPtEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
487  if(inEmcal){
488  fHistos->FillTH1(Form("hPtEMCALEtaCentNewBinning%s", eventclass), TMath::Abs(pt));
489  fHistos->FillTH1(Form("hPtEMCALEtaCentOldBinning%s", eventclass), TMath::Abs(pt));
490  }
491  for(int icut = 0; icut < 5; icut++){
492  if(TMath::Abs(pt) > static_cast<double>(ptmin[icut])){
493  fHistos->FillTH1(Form("hEtaLabDistCutPt%d%s", ptmin[icut], eventclass), etalab);
494  fHistos->FillTH1(Form("hEtaCentDistCutPt%d%s", ptmin[icut], eventclass), etacent);
495  if(inEmcal){
496  fHistos->FillTH1(Form("hEtaLabDistCutEMCALPt%d%s", ptmin[icut], eventclass), etalab);
497  fHistos->FillTH1(Form("hEtaCentDistCutEMCALPt%d%s", ptmin[icut], eventclass), etacent);
498  }
499  }
500  }
501  }
502 }
503 
509 {
510  // Called when file changes.
511 
512  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
513  if (!tree) {
514  AliError(Form("%s - UserNotify: No current tree!",GetName()));
515  return kFALSE;
516  }
517 
518  Float_t xsection = 0;
519  Float_t trials = 0;
520  Int_t pthardbin = 0;
521 
522  TFile *curfile = tree->GetCurrentFile();
523  if (!curfile) {
524  AliError(Form("%s - UserNotify: No current file!",GetName()));
525  return kFALSE;
526  }
527 
528  TChain *chain = dynamic_cast<TChain*>(tree);
529  if (chain) tree = chain->GetTree();
530 
531  Int_t nevents = tree->GetEntriesFast();
532 
533  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
534 
535  fHistos->FillTH1("hNtrials", 1., trials);
536  fHistos->FillProfile("hCrossSection", 1., xsection);
537  //fHistos->FillTH1(pthardbin, nevents);
538 
539  return kTRUE;
540 }
541 
554 Bool_t AliAnalysisTaskChargedParticlesRefMC::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard) const {
555 
556  TString file(currFile);
557  fXsec = 0;
558  fTrials = 1;
559 
560  if (file.Contains(".zip#")) {
561  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
562  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
563  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
564  file.Replace(pos+1,pos2-pos1,"");
565  } else {
566  // not an archive take the basename....
567  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
568  }
569  AliDebug(1,Form("File name: %s",file.Data()));
570 
571  // Get the pt hard bin
572  TString strPthard(file);
573 
574  strPthard.Remove(strPthard.Last('/'));
575  strPthard.Remove(strPthard.Last('/'));
576  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
577  strPthard.Remove(0,strPthard.Last('/')+1);
578  if (strPthard.IsDec())
579  pthard = strPthard.Atoi();
580  else
581  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
582 
583  // 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
584  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
585 
586  if (!fxsec) {
587  // next trial fetch the histgram file
588  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
589  if (!fxsec) {
590  // not a severe condition but inciate that we have no information
591  return kFALSE;
592  } else {
593  // find the tlist we want to be independtent of the name so use the Tkey
594  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
595  if (!key) {
596  fxsec->Close();
597  return kFALSE;
598  }
599  TList *list = dynamic_cast<TList*>(key->ReadObj());
600  if (!list) {
601  fxsec->Close();
602  return kFALSE;
603  }
604  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
605  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
606  fxsec->Close();
607  }
608  } else { // no tree pyxsec.root
609  TTree *xtree = (TTree*)fxsec->Get("Xsection");
610  if (!xtree) {
611  fxsec->Close();
612  return kFALSE;
613  }
614  UInt_t ntrials = 0;
615  Double_t xsection = 0;
616  xtree->SetBranchAddress("xsection",&xsection);
617  xtree->SetBranchAddress("ntrials",&ntrials);
618  xtree->GetEntry(0);
619  fTrials = ntrials;
620  fXsec = xsection;
621  fxsec->Close();
622  }
623  return kTRUE;
624 }
625 
631  AliGenPythiaEventHeader *pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
632  if (!pythiaHeader) {
633  // Check if AOD
634  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
635  if (aodMCH) {
636  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
637  pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
638  if (pythiaHeader) break;
639  }
640  }
641  }
642  return pythiaHeader;
643 }
644 
650  std::vector<double> mybinning;
651  std::map<double,double> definitions;
652  definitions.insert(std::pair<double,double>(2.5, 0.1));
653  definitions.insert(std::pair<double,double>(7., 0.25));
654  definitions.insert(std::pair<double,double>(15., 0.5));
655  definitions.insert(std::pair<double,double>(25., 1.));
656  definitions.insert(std::pair<double,double>(40., 2.5));
657  definitions.insert(std::pair<double,double>(50., 5.));
658  definitions.insert(std::pair<double,double>(100., 10.));
659  double currentval = 0;
660  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
661  double limit = id->first, binwidth = id->second;
662  while(currentval < limit){
663  currentval += binwidth;
664  mybinning.push_back(currentval);
665  }
666  }
667  binning.Set(mybinning.size());
668  int ib = 0;
669  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
670  binning[ib++] = *it;
671 }
672 
678  std::vector<double> mybinning;
679  std::map<double,double> definitions;
680  definitions.insert(std::pair<double, double>(1, 0.05));
681  definitions.insert(std::pair<double, double>(2, 0.1));
682  definitions.insert(std::pair<double, double>(4, 0.2));
683  definitions.insert(std::pair<double, double>(7, 0.5));
684  definitions.insert(std::pair<double, double>(16, 1));
685  definitions.insert(std::pair<double, double>(36, 2));
686  definitions.insert(std::pair<double, double>(40, 4));
687  definitions.insert(std::pair<double, double>(50, 5));
688  definitions.insert(std::pair<double, double>(100, 10));
689  definitions.insert(std::pair<double, double>(200, 20));
690  double currentval = 0.;
691  mybinning.push_back(currentval);
692  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
693  double limit = id->first, binwidth = id->second;
694  while(currentval < limit){
695  currentval += binwidth;
696  mybinning.push_back(currentval);
697  }
698  }
699  binning.Set(mybinning.size());
700  int ib = 0;
701  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
702  binning[ib++] = *it;
703 }
704 
711  return fTrackCuts->AcceptTrack(track);
712 }
713 
720  if(!track->TestFilterBit(AliAODTrack::kTrkGlobal)) return false;
721  if(track->GetTPCNCrossedRows() < 120) return false;
722  return true;
723 }
724 
730 TString AliAnalysisTaskChargedParticlesRefMC::GetFiredTriggerClasses(const TClonesArray* triggerpatches) {
731  TString triggerstring = "";
732  Int_t nEJ1 = 0, nEJ2 = 0, nEG1 = 0, nEG2 = 0;
733  double minADC_EJ1 = 260.,
734  minADC_EJ2 = 127.,
735  minADC_EG1 = 140.,
736  minADC_EG2 = 89.;
737  for(TIter patchIter = TIter(triggerpatches).Begin(); patchIter != TIter::End(); ++patchIter){
738  AliEmcalTriggerPatchInfo *patch = dynamic_cast<AliEmcalTriggerPatchInfo *>(*patchIter);
739  if(!patch->IsOfflineSimple()) continue;
740  if(patch->IsJetHighSimple() && patch->GetADCOfflineAmp() > minADC_EJ1) nEJ1++;
741  if(patch->IsJetLowSimple() && patch->GetADCOfflineAmp() > minADC_EJ2) nEJ2++;
742  if(patch->IsGammaHighSimple() && patch->GetADCOfflineAmp() > minADC_EG1) nEG1++;
743  if(patch->IsGammaLowSimple() && patch->GetADCOfflineAmp() > minADC_EG2) nEG2++;
744  }
745  if(nEJ1) triggerstring += "EJ1";
746  if(nEJ2){
747  if(triggerstring.Length()) triggerstring += ",";
748  triggerstring += "EJ2";
749  }
750  if(nEG1){
751  if(triggerstring.Length()) triggerstring += ",";
752  triggerstring += "EG1";
753  }
754  if(nEG2){
755  if(triggerstring.Length()) triggerstring += ",";
756  triggerstring += "EG2";
757  }
758  return triggerstring;
759 }
760 
769 Bool_t AliAnalysisTaskChargedParticlesRefMC::IsPhysicalPrimary(const AliVParticle* const part, AliMCEvent* const mcevent) {
770  Bool_t physprim = false;
771  const AliAODMCParticle *aodmc = dynamic_cast<const AliAODMCParticle *>(part);
772  if(aodmc){
773  physprim = aodmc->IsPhysicalPrimary();
774  } else {
775  physprim = mcevent->IsPhysicalPrimary(part->GetLabel());
776  }
777  return physprim;
778 }
779 
780 } /* namespace EMCalTriggerPtAnalysis */
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
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.
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
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)
void FillTrackHistos(const char *eventclass, Double_t pt, Double_t eta, Double_t etacent, Double_t phi, Bool_t etacut, Bool_t inEmcal)
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.
void FillTH1(const char *hname, double x, double weight=1.)
void FillProfile(const char *name, double x, double y, double weight=1.)