AliPhysics  7273240 (7273240)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskChargedParticlesRef.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 <algorithm>
16 #include <array>
17 #include <functional>
18 #include <iostream>
19 #include <map>
20 
21 #include <TMath.h>
22 #include <THistManager.h>
23 #include <TLinearBinning.h>
24 
25 #include "AliAnalysisUtils.h"
26 #include "AliAODInputHandler.h"
27 #include "AliAODTrack.h"
29 #include "AliEMCALGeometry.h"
30 #include "AliEMCALRecoUtils.h"
31 #include "AliEMCALTriggerPatchInfo.h"
32 #include "AliEmcalTrackSelection.h"
33 #include "AliESDtrackCuts.h"
34 #include "AliESDEvent.h"
35 #include "AliInputEventHandler.h"
36 #include "AliPIDResponse.h"
37 #include "AliTOFPIDResponse.h"
38 #include "AliVVertex.h"
39 
42 
46 
47 namespace EMCalTriggerPtAnalysis {
48 
52 AliAnalysisTaskChargedParticlesRef::AliAnalysisTaskChargedParticlesRef() :
54  fTrackCuts(nullptr),
55  fYshift(0.465),
56  fEtaSign(1),
57  fEtaLabCut(-0.5, 0.5),
58  fEtaCmsCut(-0.13, 0.13),
59  fPhiCut(0., TMath::TwoPi()),
60  fKineCorrelation(false),
61  fStudyPID(false)
62 {
63 }
64 
71  fTrackCuts(nullptr),
72  fYshift(0.465),
73  fEtaSign(1),
74  fEtaLabCut(-0.5, 0.5),
75  fEtaCmsCut(-0.13, 0.13),
76  fPhiCut(0., TMath::TwoPi()),
77  fKineCorrelation(false),
78  fStudyPID(false)
79 {
80  SetNeedEmcalGeom(true);
81  SetCaloTriggerPatchInfoName("EmcalTriggers");
82 }
83 
88  //if(fTrackCuts) delete fTrackCuts;
89 }
90 
92  if(!fTrackCuts) InitializeTrackCuts("standard", fInputHandler->IsA() == AliAODInputHandler::Class());
93 }
94 
99 
100  NewPtBinning newbinning;
101 
102  std::array<Double_t, 5> ptcuts = {1., 2., 5., 10., 20.};
103  // Binning for the PID histogram
104  const int kdimPID = 3;
105  const int knbinsPID[kdimPID] = {1000, 200, 300};
106  const double kminPID[kdimPID] = {-100., 0., 0.}, kmaxPID[kdimPID] = {100., 200., 1.5};
107  for(auto trg : GetSupportedTriggers()){
108  fHistos->CreateTH1(Form("hEventCount%s", trg.Data()), Form("Event Counter for trigger class %s", trg.Data()), 1, 0.5, 1.5);
109  fHistos->CreateTH1(Form("hVertexBefore%s", trg.Data()), Form("Vertex distribution before z-cut for trigger class %s", trg.Data()), 500, -50, 50);
110  fHistos->CreateTH1(Form("hVertexAfter%s", trg.Data()), Form("Vertex distribution after z-cut for trigger class %s", trg.Data()), 100, -10, 10);
111  fHistos->CreateTH1(Form("hPtEtaAll%s", trg.Data()), Form("Charged particle pt distribution all eta new binning trigger %s", trg.Data()), newbinning);
112  fHistos->CreateTH1(Form("hPtEtaCent%s", trg.Data()), Form("Charged particle pt distribution central eta new binning trigger %s", trg.Data()), newbinning);
113  fHistos->CreateTH1(Form("hPtEMCALEtaAll%s", trg.Data()), Form("Charged particle in EMCAL pt distribution all eta new binning trigger %s", trg.Data()), newbinning);
114  fHistos->CreateTH1(Form("hPtEMCALEtaCent%s", trg.Data()), Form("Charged particle in EMCAL pt distribution central eta new binning trigger %s", trg.Data()), newbinning);
115  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaAll%s", trg.Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution all eta new binning trigger %s", trg.Data()), newbinning);
116  fHistos->CreateTH1(Form("hPtEMCALNoTRDEtaCent%s", trg.Data()), Form("Charged particle in EMCAL (no TRD in front) pt distribution central eta new binning trigger %s", trg.Data()), newbinning);
117  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaAll%s", trg.Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution all eta new binning trigger %s", trg.Data()), newbinning);
118  fHistos->CreateTH1(Form("hPtEMCALWithTRDEtaCent%s", trg.Data()), Form("Charged particle in EMCAL (with TRD in front) pt distribution central eta new binning trigger %s", trg.Data()), newbinning);
119  if(fKineCorrelation){
120  fHistos->CreateTH3(Form("hPtEtaPhiAll%s", trg.Data()), Form("p_{t}-#eta-#phi distribution of all accepted tracks for trigger %s; p_{t} (GeV/c); #eta; #phi", trg.Data()), newbinning, TLinearBinning(64, -0.8, 0.8), TLinearBinning(100, 0., 2*TMath::Pi()));
121  fHistos->CreateTH3(Form("hPtEtaPhiEMCALAll%s", trg.Data()), Form("p_{t}-#eta-#phi distribution of all accepted tracks pointing to the EMCAL for trigger %s; p_{t} (GeV/c); #eta; #phi", trg.Data()), newbinning, TLinearBinning(64, -0.8, 0.8), TLinearBinning(100, 0., 2*TMath::Pi()));
122  }
123  if(fStudyPID){
124  fHistos->CreateTH2(Form("hTPCdEdxEMCAL%s", trg.Data()), Form("TPC dE/dx of charged particles in the EMCAL region for trigger %s", trg.Data()), 400, -20., 20., 200, 0., 200.);
125  fHistos->CreateTH2(Form("hTOFBetaEMCAL%s", trg.Data()), Form("TOF beta of charged particles in the EMCAL region for trigger %s", trg.Data()), 400, -20., 20., 150, 0., 1.5);
126  fHistos->CreateTHnSparse(Form("hPIDcorrEMCAL%s", trg.Data()), Form("Correlation of PID observables for Trigger %s", trg.Data()), kdimPID, knbinsPID, kminPID, kmaxPID);
127  }
128  for(auto ptcut : ptcuts) {
130  Form("hEtaLabDistAllPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
131  Form("Eta (lab) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
132  100,
133  -1.,
134  1.
135  );
137  Form("hEtaLabDistCutPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
138  Form("Eta (lab) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
139  100,
140  -1.,
141  1.
142  );
144  Form("hEtaCentDistAllPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
145  Form("Eta (cent) distribution without etacut for tracks with Pt above %.1f GeV/c trigger %s",
146  ptcut, trg.Data()),
147  160,
148  -1.3,
149  1.3
150  );
152  Form("hEtaCentDistCutPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
153  Form("Eta (cent) distribution with etacut for tracks with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
154  160,
155  -1.3,
156  1.3
157  );
159  Form("hEtaLabDistAllEMCALPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
160  Form("Eta (lab) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
161  100,
162  -1.,
163  1.
164  );
166  Form("hEtaLabDistCutEMCALPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
167  Form("Eta (lab) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
168  100,
169  -1.,
170  1.
171  );
173  Form("hEtaCentDistAllEMCALPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
174  Form("Eta (cent) distribution without etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
175  160,
176  -1.3,
177  1.3
178  );
180  Form("hEtaCentDistCutEMCALPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
181  Form("Eta (cent) distribution with etacut for tracks in EMCAL with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
182  160,
183  -1.3,
184  1.3
185  );
187  Form("hPhiDistAllPt%d%s", static_cast<Int_t>(ptcut), trg.Data()),
188  Form("#phi distribution of particles with Pt above %.1f GeV/c trigger %s", ptcut, trg.Data()),
189  300,
190  0.,
191  2*TMath::Pi()
192  );
193  }
194  }
195 }
196 
206  Bool_t hasPIDresponse = fInputHandler->GetPIDResponse() != nullptr;
207  if(fStudyPID && !hasPIDresponse) AliErrorStream() << "PID requested but PID response not available" << std::endl;
208 
209 
210  // Loop over tracks, fill select particles
211  // Histograms
212  // - Full eta_{lab} (-0.8, 0.8), new binning
213  // - Full eta_{lab} (-0.8, 0.8), old binning
214  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c without eta cut
215  // - Central eta_{cms} (-0.3, -0.2), new binning,
216  // - Central eta_{cms} (-0.8, -0.2), old binning,
217  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c
218  // - Eta distribution for tracks above 1, 2, 5, 10 GeV/c with eta cut
219  AliVTrack *checktrack(nullptr);
220  int ptmin[5] = {1,2,5,10,20}; // for eta distributions
221  Bool_t isEMCAL(kFALSE), hasTRD(kFALSE);
222  Double_t etaEMCAL(0.), phiEMCAL(0.);
223  for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
224  checktrack = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(itrk));
225  if(!checktrack) continue;
226  if(!fEtaLabCut.IsInRange(checktrack->Eta())) continue;
227  if(!fPhiCut.IsInRange(checktrack->Phi())) continue;
228  if(TMath::Abs(checktrack->Pt()) < 0.1) continue;
229  if(checktrack->IsA() == AliESDtrack::Class()){
230  AliESDtrack copytrack(*(static_cast<AliESDtrack *>(checktrack)));
231  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
232  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
233  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
234  } else {
235  AliAODTrack copytrack(*(static_cast<AliAODTrack *>(checktrack)));
236  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&copytrack);
237  etaEMCAL = copytrack.GetTrackEtaOnEMCal();
238  phiEMCAL = copytrack.GetTrackPhiOnEMCal();
239  }
240  Int_t supermoduleID = -1;
241  isEMCAL = fGeom->SuperModuleNumberFromEtaPhi(etaEMCAL, phiEMCAL, supermoduleID);
242  // Exclude supermodules 10 and 11 as they did not participate in the trigger
243  isEMCAL = isEMCAL && supermoduleID < 10;
244  hasTRD = isEMCAL && supermoduleID >= 4; // supermodules 4 - 10 have TRD in front in the 2012-2013 ALICE setup
245 
246  // Calculate eta in cms frame according
247  // EPJC74 (2014) 3054:
248  // eta_cms = - eta_lab - |yshift|
249  Double_t etacent = -1. * checktrack->Eta() - TMath::Abs(fYshift);
250  etacent *= fEtaSign;
251 
252  Bool_t etacentcut = fEtaCmsCut.IsInRange(etacent);
253 
254  if(!fTrackCuts->IsTrackAccepted(checktrack)) continue;
255 
256  for(const auto &t : fSelectedTriggers){
257  FillTrackHistos(t, checktrack->Pt(), checktrack->Eta() * fEtaSign, etacent, checktrack->Phi(), etacentcut, isEMCAL, hasTRD);
258  if(fStudyPID && hasPIDresponse)
259  if(isEMCAL) FillPIDHistos(t, *checktrack);
260  }
261  }
262  return true;
263 }
264 
266  // Apply vertex z cut
267  for(const auto &t : fSelectedTriggers){
268  Double_t weight = GetTriggerWeight(t);
269  fHistos->FillTH1(Form("hVertexBefore%s", t.Data()), fVertex[2], weight);
270  }
271 }
272 
274  for(const auto &t : fSelectedTriggers) {
275  Double_t weight = GetTriggerWeight(t);
276  // Fill Event counter and reference vertex distributions after event selection
277  fHistos->FillTH1(Form("hEventCount%s", t.Data()), 1, weight);
278  fHistos->FillTH1(Form("hVertexAfter%s", t.Data()), fVertex[2], weight);
279  }
280 }
281 
282 
294  const TString &eventclass,
295  Double_t pt,
296  Double_t etalab,
297  Double_t etacent,
298  Double_t phi,
299  Bool_t etacut,
300  Bool_t inEmcal,
301  Bool_t hasTRD
302  )
303 {
304  Double_t weight = GetTriggerWeight(eventclass);
305  AliDebugStream(1) << GetName() << ": Using weight " << weight << " for trigger " << eventclass << " in particle histograms." << std::endl;
306  fHistos->FillTH1(Form("hPtEtaAll%s", eventclass.Data()), TMath::Abs(pt), weight);
307  double kinepoint[3] = {TMath::Abs(pt), etalab, phi};
308  if(fKineCorrelation) fHistos->FillTH3(Form("hPtEtaPhiAll%s", eventclass.Data()),kinepoint , weight);
309  if(inEmcal){
310  fHistos->FillTH1(Form("hPtEMCALEtaAll%s", eventclass.Data()), TMath::Abs(pt), weight);
311  if(fKineCorrelation) fHistos->FillTH3(Form("hPtEtaPhiEMCALAll%s", eventclass.Data()), kinepoint, weight);
312  if(hasTRD){
313  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaAll%s", eventclass.Data()), TMath::Abs(pt), weight);
314  } else {
315  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaAll%s", eventclass.Data()), TMath::Abs(pt), weight);
316  }
317  }
318 
319  std::array<int, 5> ptmin = {1,2,5,10,20}; // for eta distributions
320  for(auto ptmincut : ptmin){
321  if(TMath::Abs(pt) > static_cast<double>(ptmincut)){
322  fHistos->FillTH1(Form("hPhiDistAllPt%d%s", ptmincut, eventclass.Data()), phi, weight);
323  fHistos->FillTH1(Form("hEtaLabDistAllPt%d%s", ptmincut, eventclass.Data()), etalab, weight);
324  fHistos->FillTH1(Form("hEtaCentDistAllPt%d%s", ptmincut, eventclass.Data()), etacent, weight);
325  if(inEmcal){
326  fHistos->FillTH1(Form("hEtaLabDistAllEMCALPt%d%s", ptmincut, eventclass.Data()), etalab, weight);
327  fHistos->FillTH1(Form("hEtaCentDistAllEMCALPt%d%s", ptmincut, eventclass.Data()), etacent, weight);
328  }
329  }
330  }
331 
332  if(etacut){
333  fHistos->FillTH1(Form("hPtEtaCent%s", eventclass.Data()), TMath::Abs(pt), weight);
334  if(inEmcal){
335  fHistos->FillTH1(Form("hPtEMCALEtaCent%s", eventclass.Data()), TMath::Abs(pt), weight);
336  if(hasTRD){
337  fHistos->FillTH1(Form("hPtEMCALWithTRDEtaCent%s", eventclass.Data()), TMath::Abs(pt), weight);
338  } else {
339  fHistos->FillTH1(Form("hPtEMCALNoTRDEtaCent%s", eventclass.Data()), TMath::Abs(pt), weight);
340  }
341  }
342  for(auto ptmincut : ptmin){
343  if(TMath::Abs(pt) > static_cast<double>(ptmincut)){
344  fHistos->FillTH1(Form("hEtaLabDistCutPt%d%s", ptmincut, eventclass.Data()), etalab, weight);
345  fHistos->FillTH1(Form("hEtaCentDistCutPt%d%s", ptmincut, eventclass.Data()), etacent, weight);
346  if(inEmcal){
347  fHistos->FillTH1(Form("hEtaLabDistCutEMCALPt%d%s", ptmincut, eventclass.Data()), etalab, weight);
348  fHistos->FillTH1(Form("hEtaCentDistCutEMCALPt%d%s", ptmincut, eventclass.Data()), etacent, weight);
349  }
350  }
351  }
352  }
353 }
354 
356  const TString &eventclass,
357  const AliVTrack &trk
358 ) {
359  Double_t weight = GetTriggerWeight(eventclass);
360  AliDebugStream(1) << GetName() << ": Using weight " << weight << " for trigger " << eventclass << " in PID histograms." << std::endl;
361  AliPIDResponse *pid = fInputHandler->GetPIDResponse();
362  if(TMath::Abs(trk.Eta()) > 0.5) return;
363  if(!((trk.GetStatus() & AliVTrack::kTOFout) && (trk.GetStatus() & AliVTrack::kTIME))) return;
364 
365  double poverz = TMath::Abs(trk.P())/static_cast<double>(trk.Charge());
366  fHistos->FillTH2(Form("hTPCdEdxEMCAL%s", eventclass.Data()), poverz, trk.GetTPCsignal(), weight);
367  // correct for units - TOF in ps, track length in cm
368  Double_t trtime = (trk.GetTOFsignal() - pid->GetTOFResponse().GetTimeZero()) * 1e-12;
369  Double_t v = trk.GetIntegratedLength()/(100. * trtime);
370  Double_t beta = v / TMath::C();
371  fHistos->FillTH2(Form("hTOFBetaEMCAL%s", eventclass.Data()), poverz, beta, weight);
372  double datapoint[3] = {poverz, trk.GetTPCsignal(), beta};
373  fHistos->FillTHnSparse(Form("hPIDcorrEMCAL%s", eventclass.Data()), datapoint, weight);
374 }
375 
376 
384 }
385 
392 {
393  this->SetMinimum(0.);
394  this->AddStep(1, 0.05);
395  this->AddStep(2, 0.1);
396  this->AddStep(4, 0.2);
397  this->AddStep(7, 0.5);
398  this->AddStep(16, 1);
399  this->AddStep(36, 2);
400  this->AddStep(40, 4);
401  this->AddStep(50, 5);
402  this->AddStep(100, 10);
403  this->AddStep(200, 20);
404 }
405 
406 } /* namespace EMCalTriggerPtAnalysis */
std::vector< TString > fSelectedTriggers
! Triggers selected for given event
double Double_t
Definition: External.C:58
Double_t fEtaSign
Sign of the eta distribution (swaps when beam directions swap): p-Pb: +1, Pb-p: -1.
AliEmcalTrackSelection * fTrackCuts
Standard track selection.
Class creating a linear binning, used in the histogram manager.
static AliEmcalTrackSelection * TrackCutsFactory(TString name, Bool_t isAOD)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
void AddStep(Double_t max, Double_t binwidth)
void SetCaloTriggerPatchInfoName(const char *n)
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
int Int_t
Definition: External.C:63
AliEMCALGeometry * fGeom
!emcal geometry
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Helper class creating user defined custom binning.
const Double_t ptmin
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Declaration of class AliEMCalTriggerExtraCuts.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t fVertex[3]
!event vertex
void FillPIDHistos(const TString &eventclass, const AliVTrack &track)
void SetNeedEmcalGeom(Bool_t n)
bool Bool_t
Definition: External.C:53
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
TH3 * CreateTH3(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, int nbinsz, double zmin, double zmax, Option_t *opt="")
virtual bool IsTrackAccepted(AliVTrack *const trk)=0
void FillTrackHistos(const TString &eventclass, Double_t pt, Double_t eta, Double_t etacent, Double_t phi, Bool_t etacut, Bool_t inEmcal, Bool_t hasTRD)
void SetMinimum(Double_t min)