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