AliPhysics  v5-06-21-01 (1eac791)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliAnalysisTaskPtEMCalTrigger.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2014, 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 
16 /*
17  * Analysis task of the pt analysis on EMCal-triggered events
18  *
19  * Author: Markus Fasel
20  */
21 
22 #include <map>
23 #include <cstring>
24 #include <iostream>
25 #include <memory>
26 #include <vector>
27 #include <string>
28 #include <sstream>
29 
30 #include <TClonesArray.h>
31 #include <TDirectory.h>
32 #include <TH1.h>
33 #include <THashList.h>
34 #include <TList.h>
35 #include <TKey.h>
36 #include <TMath.h>
37 #include <TObjArray.h>
38 #include <TObjString.h>
39 #include <TString.h>
40 #include <TVector2.h>
41 
42 #include "AliAODEvent.h"
43 #include "AliESDEvent.h"
44 #include "AliInputEventHandler.h"
45 #include "AliMCEvent.h"
46 #include "AliParticleContainer.h"
47 #include "AliVCluster.h"
48 #include "AliVParticle.h"
49 #include "AliVTrack.h"
50 #include "AliVVertex.h"
51 
52 #include "AliClusterContainer.h"
53 #include "AliEmcalJet.h"
54 #include "AliEmcalTriggerPatchInfo.h"
55 #include "AliEmcalPhysicsSelection.h"
56 #include "AliAnalysisTaskPtEMCalTrigger.h"
57 #include "AliEMCalHistoContainer.h"
58 #include "AliEMCalPtTaskVTrackSelection.h"
59 #include "AliEMCalPtTaskTrackSelectionESD.h"
60 #include "AliEMCalPtTaskTrackSelectionAOD.h"
61 #include "AliJetContainer.h"
62 #include "AliParticleContainer.h"
63 #include "AliPicoTrack.h"
64 
65 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)
66 
67 namespace EMCalTriggerPtAnalysis {
68 
69  /*
70  * constants
71  */
72  const Int_t AliAnalysisTaskPtEMCalTrigger::kNJetRadii = 4;
73  const Double_t jetRadVals[AliAnalysisTaskPtEMCalTrigger::kNJetRadii] = {0.2, 0.3, 0.4, 0.5};
74  const Double_t *AliAnalysisTaskPtEMCalTrigger::kJetRadii = jetRadVals;
75 
76  //______________________________________________________________________________
77  AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
78  AliAnalysisTaskEmcalJet(),
79  fHistos(NULL),
80  fListTrackCuts(NULL),
81  fEtaRange(),
82  fPtRange(),
83  fEnergyRange(),
84  fVertexRange(),
85  fJetContainersMC(),
86  fJetContainersData(),
87  fSelectAllTracks(kFALSE),
88  fSwapEta(kFALSE),
89  fUseTriggersFromTriggerMaker(kFALSE)
90  {
91  /*
92  * Dummy constructor, initialising the values with default (NULL) values
93  */
94  }
95 
96  //______________________________________________________________________________
97  AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
98  AliAnalysisTaskEmcalJet(name, kTRUE),
99  fHistos(NULL),
100  fListTrackCuts(NULL),
101  fEtaRange(),
102  fPtRange(),
103  fEnergyRange(),
104  fVertexRange(),
105  fJetContainersMC(),
106  fJetContainersData(),
107  fSelectAllTracks(kFALSE),
108  fSwapEta(kFALSE),
109  fUseTriggersFromTriggerMaker(kFALSE)
110  {
111  /*
112  * Main constructor, setting default values for eta and zvertex cut
113  */
114 
115  fListTrackCuts = new TList;
116  fListTrackCuts->SetOwner(false);
117 
118  // Set default cuts
119  fEtaRange.SetLimits(-0.8, 0.8);
120  fPtRange.SetLimits(0.15, 100.);
121  fEnergyRange.SetLimits(0., 1000.);
122  fVertexRange.SetLimits(-10., 10.);
123  SetMakeGeneralHistograms(kTRUE);
124  }
125 
126  //______________________________________________________________________________
127  AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
128  /*
129  * Destructor, deleting output
130  */
131  //if(fTrackSelection) delete fTrackSelection;
132  if(fHistos) delete fHistos;
133  if(fListTrackCuts) delete fListTrackCuts;
134  }
135 
136  //______________________________________________________________________________
137  void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
138  /*
139  * Create the list of output objects and define the histograms.
140  * Also adding the track cuts to the list of histograms.
141  */
142  AliAnalysisTaskEmcal::UserCreateOutputObjects();
143  SetCaloTriggerPatchInfoName("EmcalTriggers");
144  fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
145  fHistos->ReleaseOwner();
146 
147  if(fJetContainersMC.GetEntries()){
148  AliDebug(1,"Jet containers for MC truth:");
149  TObjString *contname(NULL);
150  TIter contMCIter(&fJetContainersMC);
151  while((contname = dynamic_cast<TObjString *>(contMCIter.Next())))
152  AliDebug(1, Form("Next container: %s", contname->String().Data()));
153  }
154  if(fJetContainersData.GetEntries()){
155  AliDebug(1, "Jet containers for Data:");
156  TObjString *contname(NULL);
157  TIter contDataIter(&fJetContainersData);
158  while((contname = dynamic_cast<TObjString *>(contDataIter.Next())))
159  AliDebug(1, Form("Next container: %s", contname->String().Data()));
160  }
161  if(fJetCollArray.GetEntries()){
162  AliDebug(1, "Jet containers attached to this task:");
163  AliJetContainer *cont(NULL);
164  TIter contIter(&fJetCollArray);
165  while((cont = dynamic_cast<AliJetContainer *>(contIter.Next())))
166  AliDebug(1, Form("Container: %s", cont->GetName()));
167  }
168 
169  std::map<std::string, std::string> triggerCombinations;
170  const char *triggernames[12] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal", "EMCHighBoth", "EMCHighGammaOnly", "EMCHighJetOnly", "EMCLowBoth", "EMCLowGammaOnly", "EMCLowJetOnly"},
171  *bitnames[4] = {"CINT7", "EMC7", "kEMCEGA", "kEMCEJE"};
172  // Define axes for the trigger correlation histogram
173  const TAxis *triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis *) * 5);
174  const TAxis *bitaxes[4]; memset(bitaxes, 0, sizeof(TAxis *) * 4);
175  const char *binlabels[2] = {"OFF", "ON"};
176  TAxis mytrgaxis[5], mybitaxis[4];
177  for(int itrg = 0; itrg < 5; ++itrg){
178  DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
179  triggeraxis[itrg] = mytrgaxis+itrg;
180  if(itrg < 4){
181  DefineAxis(mybitaxis[itrg], bitnames[itrg], bitnames[itrg], 2, -0.5, 1.5, binlabels);
182  bitaxes[itrg] = mybitaxis+itrg;
183  }
184  }
185  // Define names and titles for different triggers in the histogram container
186  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
187  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
188  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
189  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "gamma-triggered events (high threshold)"));
190  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "gamma-triggered events (low threshold)"));
191  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events"));
192  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[6], "jet and gamma triggered events (high threshold)"));
193  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[7], "exclusively gamma-triggered events (high threshold)"));
194  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[8], "exclusively jet-triggered events (high threshold)"));
195  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[9], "jet and gamma triggered events (low threshold)"));
196  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[10], "exclusively gamma-triggered events (low threshold)"));
197  triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[11], "exclusively-triggered events (low threshold)"));
198  // Define axes for the pt histogram
199  // Dimensions:
200  // 1. pt
201  // 2. eta
202  // 3. phi
203  // 4. vertex
204  // 5. pileup (0 = all events, 1 = after pileup rejection)
205  // 6. track cuts (0 = no cuts; 1 = after std cuts)
206  TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
207  pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
208  CreateDefaultPtBinning(ptbinning);
209  CreateDefaultZVertexBinning(zvertexBinning);
210  CreateDefaultEtaBinning(etabinning);
211  TAxis htrackaxes[7];
212  DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
213  DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
214  DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
215  DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
216  DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
217  DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
218  DefineAxis(htrackaxes[6], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
219  const TAxis *trackaxes[7];
220  for(int iaxis = 0; iaxis < 7; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
221  TAxis hclusteraxes[4];
222  DefineAxis(hclusteraxes[0], "energy", "E (GeV)", ptbinning);
223  DefineAxis(hclusteraxes[1], "zvertex", "z_{V} (cm)", zvertexBinning);
224  DefineAxis(hclusteraxes[2], "pileup", "Pileup rejection", 2,1 -0.5, 1.5);
225  DefineAxis(hclusteraxes[3], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
226  const TAxis *clusteraxes[4];
227  for(int iaxis = 0; iaxis < 4; ++iaxis) clusteraxes[iaxis] = hclusteraxes + iaxis;
228  TAxis hpatchenergyaxes[5];
229  DefineAxis(hpatchenergyaxes[0], "energy", "Patch energy (GeV)", 100, 0., 100);
230  DefineAxis(hpatchenergyaxes[1], "eta", "#eta", etabinning);
231  DefineAxis(hpatchenergyaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
232  DefineAxis(hpatchenergyaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
233  DefineAxis(hpatchenergyaxes[4], "emcalgood", "EMCAL good event", 2, -0.5, 1.5);
234  const TAxis *patchenergyaxes[5];
235  for(int iaxis = 0; iaxis < 5; ++iaxis) patchenergyaxes[iaxis] = hpatchenergyaxes + iaxis;
236  TAxis hpatchampaxes[5];
237  DefineAxis(hpatchampaxes[0], "amplitude", "Patch energy (GeV)", 10000, 0., 10000.);
238  DefineAxis(hpatchampaxes[1], "eta", "#eta", etabinning);
239  DefineAxis(hpatchampaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
240  DefineAxis(hpatchampaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
241  DefineAxis(hpatchampaxes[4], "emcalgood", "EMCAL good event", 2, -0.5, 1.5);
242  const TAxis *patchampaxes[5];
243  for(int iaxis = 0; iaxis < 5; ++iaxis) patchampaxes[iaxis] = hpatchampaxes + iaxis;
244  try{
245  std::string patchnames[] = {"Level0", "JetHigh", "JetLow", "GammaHigh", "GammaLow"};
246  for(std::string * triggerpatch = patchnames; triggerpatch < patchnames + sizeof(patchnames)/sizeof(std::string); ++triggerpatch){
247  fHistos->CreateTHnSparse(Form("Energy%s", triggerpatch->c_str()), Form("Patch energy for %s trigger patches", triggerpatch->c_str()), 5, patchenergyaxes, "s");
248  fHistos->CreateTHnSparse(Form("EnergyRough%s", triggerpatch->c_str()), Form("Rough patch energy for %s trigger patches", triggerpatch->c_str()), 5, patchenergyaxes, "s");
249  fHistos->CreateTHnSparse(Form("Amplitude%s", triggerpatch->c_str()), Form("Patch amplitude for %s trigger patches", triggerpatch->c_str()), 5, patchampaxes, "s");
250  }
251 
252  // Create histogram for MC-truth
253  fHistos->CreateTHnSparse("hMCtrueParticles", "Particle-based histogram for MC-true particles", 5, trackaxes, "s");
254  if(fJetCollArray.GetEntries()){
255  for(int irad = 0; irad < kNJetRadii; irad++){
256  fHistos->CreateTHnSparse(Form("hMCtrueParticlesRad%02d", int(kJetRadii[irad]*10)),
257  Form("Particle-based histogram for MC-true particles in Jets with radius %.1f", kJetRadii[irad]*10), 5, trackaxes, "s");
258  }
259  // histogram for isolated particles
260  fHistos->CreateTHnSparse("hMCtrueParticlesIsolated", "Particle-based histogram for isolated MC-true particles", 5, trackaxes, "s");
261  }
262  for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
263  const std::string name = it->first, &title = it->second;
264  // Create event-based histogram
265  fHistos->CreateTH2(Form("hEventHist%s", name.c_str()), Form("Event-based data for %s events; pileup rejection; z_{V} (cm)", title.c_str()), pileupaxis, zvertexBinning);
266  // Create track-based histogram
267  fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes, "s");
268  fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events for tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
269  fHistos->CreateTHnSparse(Form("hMCTrackHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes, "s");
270  fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics for tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
271  // Create cluster-based histogram (Uncalibrated and calibrated clusters)
272  fHistos->CreateTHnSparse(Form("hClusterCalibHist%s", name.c_str()), Form("Calib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes, "s");
273  fHistos->CreateTHnSparse(Form("hClusterUncalibHist%s", name.c_str()), Form("Uncalib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes, "s");
274  if(fJetCollArray.GetEntries()){
275  // Create histograms for particles in different jetfinder radii, only track-based
276  for(int irad = 0; irad < kNJetRadii; irad++){
277  fHistos->CreateTHnSparse(Form("hTrackHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
278  Form("Track-based data for %s events for tracks in jets with Radius %.1f", title.c_str(), kJetRadii[irad]), 7, trackaxes, "s");
279  fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
280  Form("Track-based data for %s events for tracks matched to EMCal clusters in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
281  7, trackaxes, "s");
282  fHistos->CreateTHnSparse(Form("hMCTrackHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
283  Form("Track-based data for %s events with MC kinematics for tracks in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
284  7, trackaxes, "s");
285  fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
286  Form("Track-based data for %s events with MC kinematics for tracks matched to EMCal clusters in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
287  7, trackaxes, "s");
288  fHistos->CreateTHnSparse(Form("hClusterCalibHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
289  Form("Calib. cluster-based histogram for %s events for clusters in Jets with radius %.1f", title.c_str(), kJetRadii[irad]), 4, clusteraxes, "s");
290  }
291  }
292  // Create also histograms for isolated particles
293  fHistos->CreateTHnSparse(Form("hTrackHist%sIsolated", name.c_str()), Form("Track-based data for %s events for isolated tracks", title.c_str()), 7, trackaxes, "s");
294  fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%sIsolated", name.c_str()), Form("Track-based data for %s events for isolated tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
295  fHistos->CreateTHnSparse(Form("hMCTrackHist%sIsolated", name.c_str()), Form("Track-based data for %s events with MC kinematics for isolated tracks", title.c_str()), 7, trackaxes, "s");
296  fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%sIsolated", name.c_str()), Form("Track-based data for %s events with MC kinematics for isolated tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
297  }
298  fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
299  fHistos->CreateTHnSparse("hEventsTriggerbit", "Trigger bits for the different events", 4, bitaxes);
300  } catch (HistoContainerContentException &e){
301  std::stringstream errormessage;
302  errormessage << "Creation of histogram failed: " << e.what();
303  AliError(errormessage.str().c_str());
304  }
305  fOutput->Add(fHistos->GetListOfHistograms());
306  if(fListTrackCuts && fListTrackCuts->GetEntries()){
307  TIter cutIter(fListTrackCuts);
308  AliEMCalPtTaskVTrackSelection *cutObject(NULL);
309  while((cutObject = dynamic_cast<AliEMCalPtTaskVTrackSelection *>(cutIter()))){
310  AliESDtrackCuts *cuts = dynamic_cast<AliESDtrackCuts *>(cutObject->GetTrackCuts(0));
311  if(cuts){
312  cuts->DefineHistograms();
313  fOutput->Add(cuts);
314  }
315  }
316  }
317  PostData(1, fOutput);
318  }
319 
320  //______________________________________________________________________________
321  Bool_t AliAnalysisTaskPtEMCalTrigger::Run(){
322  /*
323  * Runs the event loop
324  *
325  * @param option: Additional options
326  */
327  // Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
328  AliDebug(1,Form("Number of calibrated clusters: %d", fCaloClusters->GetEntries()));
329  AliDebug(1,Form("Number of matched tracks: %d", fTracks->GetEntries()));
330  if(fMCEvent){
331  // Build always trigger strig from the trigger maker in case of MC
332  fUseTriggersFromTriggerMaker = kTRUE;
333  }
334 
335  Bool_t emcalGood = fInputHandler->IsEventSelected() & AliEmcalPhysicsSelection::kEmcalOk;
336 
337  // Loop over trigger patches, fill patch energy
338  AliEmcalTriggerPatchInfo *triggerpatch(NULL);
339  TIter patchIter(this->fTriggerPatchInfo);
340  while((triggerpatch = dynamic_cast<AliEmcalTriggerPatchInfo *>(patchIter()))){
341  double triggerpatchinfo[5] = {triggerpatch->GetPatchE(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
342  double triggerpatchinfoamp[5] = {static_cast<double>(triggerpatch->GetADCAmp()), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
343  double triggerpatchinfoer[5] = {triggerpatch->GetADCAmpGeVRough(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
344  if(triggerpatch->IsJetHigh()){
345  fHistos->FillTHnSparse("EnergyJetHigh", triggerpatchinfo);
346  fHistos->FillTHnSparse("AmplitudeJetHigh", triggerpatchinfoamp);
347  fHistos->FillTHnSparse("EnergyRoughJetHigh", triggerpatchinfoer);
348  }
349  if(triggerpatch->IsJetLow()){
350  fHistos->FillTHnSparse("EnergyJetLow", triggerpatchinfo);
351  fHistos->FillTHnSparse("AmplitudeJetLow", triggerpatchinfoamp);
352  fHistos->FillTHnSparse("EnergyRoughJetLow", triggerpatchinfoer);
353  }
354  if(triggerpatch->IsGammaHigh()){
355  fHistos->FillTHnSparse("EnergyGammaHigh", triggerpatchinfo);
356  fHistos->FillTHnSparse("AmplitudeGammaHigh", triggerpatchinfoamp);
357  fHistos->FillTHnSparse("EnergyRoughGammaHigh", triggerpatchinfoer);
358  }
359  if(triggerpatch->IsGammaLow()){
360  fHistos->FillTHnSparse("EnergyGammaLow", triggerpatchinfo);
361  fHistos->FillTHnSparse("AmplitudeGammaLow", triggerpatchinfoamp);
362  fHistos->FillTHnSparse("EnergyRoughGammaLow", triggerpatchinfoer);
363  }
364  if(triggerpatch->IsLevel0()){
365  fHistos->FillTHnSparse("EnergyLevel0", triggerpatchinfo);
366  fHistos->FillTHnSparse("AmplitudeLevel0", triggerpatchinfoamp);
367  fHistos->FillTHnSparse("EnergyRoughLevel0", triggerpatchinfoer);
368  }
369  }
370 
371  const AliVVertex *vtxTracks = fInputEvent->GetPrimaryVertex(),
372  *vtxSPD = GetSPDVertex();
373  if(!(vtxTracks && vtxSPD)) return false;
374  if(!fVertexRange.IsInRange(vtxTracks->GetZ())) return false;
375  if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return false;
376 
377  double triggers[5]; memset(triggers, 0, sizeof(double) * 5);
378  double triggerbits[4]; memset(triggerbits, 0, sizeof(double) * 4);
379  if(fInputHandler->IsEventSelected() & AliVEvent::kINT7){
380  triggers[0] = 1.;
381  triggerbits[0] = 1.;
382  }
383 
384  // check triggerbits
385  if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
386  triggerbits[1] = 1.;
387  }
388  if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA){
389  triggerbits[2] = 1.;
390  }
391  if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE){
392  triggerbits[3] = 1.;
393  }
394  try{
395  fHistos->FillTHnSparse("hEventsTriggerbit", triggerbits);
396  } catch(HistoContainerContentException &e) {
397  std::stringstream errormessage;
398  errormessage << "Filling of histogram failed: " << e.what();
399  AliError(errormessage.str().c_str());
400  }
401 
402  std::vector<std::string> triggerstrings;
403  // EMCal-triggered event, distinguish types
404  TString trgstr(fUseTriggersFromTriggerMaker ? BuildTriggerString() : fInputEvent->GetFiredTriggerClasses());
405  AliDebug(1, Form("Triggerstring: %s\n", trgstr.Data()));
406  if(trgstr.Contains("EJ1")){
407  triggerstrings.push_back("EMCJHigh");
408  triggers[1] = 1;
409  if(trgstr.Contains("EG1"))
410  triggerstrings.push_back("EMCHighBoth");
411  else
412  triggerstrings.push_back("EMCHighJetOnly");
413  }
414  if(trgstr.Contains("EJ2")){
415  triggerstrings.push_back("EMCJLow");
416  triggers[2] = 1;
417  if(trgstr.Contains("EG2"))
418  triggerstrings.push_back("EMCLowBoth");
419  else
420  triggerstrings.push_back("EMCLowJetOnly");
421  }
422  if(trgstr.Contains("EG1")){
423  triggerstrings.push_back("EMCGHigh");
424  triggers[3] = 1;
425  if(!trgstr.Contains("EJ1"))
426  triggerstrings.push_back("EMCHighGammaOnly");
427  }
428  if(trgstr.Contains("EG2")){
429  triggerstrings.push_back("EMCGLow");
430  triggers[4] = 1;
431  if(!trgstr.Contains("EJ2"))
432  triggerstrings.push_back("EMCLowGammaOnly");
433  }
434 
435  try{
436  fHistos->FillTHnSparse("hEventTriggers", triggers);
437  } catch (HistoContainerContentException &e){
438  std::stringstream errormessage;
439  errormessage << "Filling of histogram failed: " << e.what();
440  AliError(errormessage.str().c_str());
441  }
442 
443  // apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
444  bool isPileupEvent = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.);
445  isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
446  double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
447  isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
448 
449  // Fill event-based histogram
450  const double &zv = vtxTracks->GetZ();
451  if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
452  if(!triggerstrings.size()) // Non-EMCal-triggered
453  FillEventHist("NoEMCal", zv, isPileupEvent);
454  else{
455  // EMCal-triggered events
456  for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
457  FillEventHist(it->c_str(), zv, isPileupEvent);
458  }
459 
460  const AliEmcalJet *foundJet(NULL);
461  char histname[1024];
462  std::vector<double> coneradii;
463  AliJetContainer *jetContMC(NULL), *jetContData(NULL);
464 
465  // Fill MC truth
466  AliVParticle *part(NULL);
467  if(fMCEvent){
468  if(fJetCollArray.GetEntries() &&
469  (jetContMC = dynamic_cast<AliJetContainer *>(fJetCollArray.FindObject((static_cast<TObjString *>(fJetContainersMC.At(0)))->String().Data())))){
470  // In case we have a jet array we loop over all MC selected particles in the particle container of the jet array
471  TIter particles(jetContMC->GetParticleContainer()->GetArray());
472  while((part = dynamic_cast<AliVParticle *>(particles()))){
473  if(part->Charge() == 0) continue;
474  if(!fEtaRange.IsInRange(part->Eta())) continue;
475  if(!fPtRange.IsInRange(part->Pt())) continue;
476  FillMCParticleHist("hMCtrueParticles", part, zv, isPileupEvent);
477 
478  /*
479  * Jet part: Find track in jet container,
480  * check according to number of particles in jet, and
481  * check for different cone radii
482  */
483  foundJet = FoundTrackInJet(part, jetContMC);
484  if(foundJet && foundJet->GetNumberOfConstituents() > 1){
485  for(int irad = 0; irad < kNJetRadii; irad++){
486  if(IsInRadius(part, foundJet, kJetRadii[irad])){
487  sprintf(histname, "hMCtrueParticlesRad%02d", int(kJetRadii[irad]*10));
488  FillMCParticleHist(histname, part, zv, isPileupEvent);
489  }
490  }
491  } else {
492  // isolated track
493  FillMCParticleHist("hMCtrueParticlesIsolated", part, zv, isPileupEvent);
494  }
495  }
496  } else {
497  // Use MC Event
498  for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
499  // Select only physical primary particles
500  part = fMCEvent->GetTrack(ipart);
501  if(part->Charge() == 0) continue;
502  if(!fEtaRange.IsInRange(part->Eta())) continue;
503  if(!fPtRange.IsInRange(part->Pt())) continue;
504  if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
505  FillMCParticleHist("hMCtrueParticles", part, zv, isPileupEvent);
506  }
507  }
508  }
509 
510  AliVTrack *track(NULL);
511  AliPicoTrack *picoTrack(NULL);
512  TObject *containerObject(NULL);
513  // Loop over all tracks (No cuts applied)
514  if(fSelectAllTracks){
515  // loop over all tracks only if requested
516  TIter allTrackIter(fTracks);
517  while((containerObject = dynamic_cast<TObject *>(allTrackIter()))){
518  if((picoTrack = dynamic_cast<AliPicoTrack *>(containerObject))){
519  track = picoTrack->GetTrack();
520  } else
521  track = dynamic_cast<AliVTrack *>(containerObject);
522  if(!IsTrueTrack(track)) continue;
523  if(!fEtaRange.IsInRange(track->Eta())) continue;
524  if(!fPtRange.IsInRange(track->Pt())) continue;
525  if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0, triggers[0]);
526  if(!triggerstrings.size()) // Non-EMCal-triggered
527  FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0, triggers[0]);
528  else {
529  // EMCal-triggered events
530  for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
531  FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0, triggers[0]);
532  }
533  }
534  }
535 
536  // Now apply track selection cuts
537  // allow for several track selections to be tested at the same time
538  // each track selection gets a different cut ID starting from 1
539  // cut ID 0 is reserved for the case of no cuts
540  // In case we have a jet container attached we check whether the track is
541  // found in a jet, and check for different cone radii around a jet
542  if(fListTrackCuts && fListTrackCuts->GetEntries()){
543  for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
544  AliEMCalPtTaskVTrackSelection *trackSelection = static_cast<AliEMCalPtTaskVTrackSelection *>(fListTrackCuts->At(icut));
545  TIter trackIter(trackSelection->GetAcceptedTracks(fTracks));
546  while((track = dynamic_cast<AliVTrack *>(trackIter()))){
547  if(fMCEvent && !IsTrueTrack(track)) continue; // Reject fake tracks in case of MC
548  if(!fEtaRange.IsInRange(track->Eta())) continue;
549  if(!fPtRange.IsInRange(track->Pt())) continue;
550  coneradii.clear();
551  if(this->fJetCollArray.GetEntries() &&
552  (jetContData = dynamic_cast<AliJetContainer *>(this->fJetCollArray.FindObject((static_cast<TObjString *>(this->fJetContainersData.At(0)))->String().Data())))){
553  foundJet = FoundTrackInJet(track, jetContData);
554  if(foundJet){
555  for(int irad = 0; irad < kNJetRadii; irad++){
556  if(IsInRadius(track, foundJet, kJetRadii[irad])) coneradii.push_back(kJetRadii[irad]);
557  }
558  }
559  }
560  if(triggers[0]){
561  FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0]);
562  if(coneradii.size()){
563  for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
564  FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
565  }
566  }
567  if(!triggerstrings.size()){ // Non-EMCal-triggered
568  FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0]);
569  for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
570  FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
571  } else {
572  // EMCal-triggered events
573  for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
574  FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0]);
575  for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
576  FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
577  }
578  }
579  }
580  }
581  }
582 
583  // Next step we loop over the (uncalibrated) emcal clusters and fill histograms with the cluster energy
584  const AliVCluster *clust(NULL);
585  for(int icl = 0; icl < fInputEvent->GetNumberOfCaloClusters(); icl++){
586  clust = fInputEvent->GetCaloCluster(icl);
587  if(!clust->IsEMCAL()) continue;
588  if(!fEnergyRange.IsInRange(clust->E())) continue;
589  if(triggers[0]) FillClusterHist("hClusterUncalibHistMinBias", clust, zv, isPileupEvent, triggers[0]);
590  if(!triggerstrings.size()){ // Non-EMCal-triggered
591  FillClusterHist("hClusterUncalibHistNoEMCal", clust, zv, isPileupEvent, triggers[0]);
592  } else{
593  for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
594  sprintf(histname, "hClusterUncalibHist%s", it->c_str());
595  FillClusterHist(histname, clust, zv, isPileupEvent, triggers[0]);
596  }
597  }
598  }
599 
600  if(fCaloClusters){
601  TIter clustIter(fCaloClusters);
602  while((clust = dynamic_cast<const AliVCluster *>(clustIter()))){
603  if(!clust->IsEMCAL()) continue;
604  if(!fEnergyRange.IsInRange(clust->E())) continue;
605  if(triggers[0]) FillClusterHist("hClusterCalibHistMinBias", clust, zv, isPileupEvent, triggers[0]);
606  if(!triggerstrings.size()) // Non-EMCal-triggered
607  FillClusterHist("hClusterCalibHistNoEMCal", clust, zv, isPileupEvent, triggers[0]);
608  else{
609  for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
610  sprintf(histname, "hClusterCalibHist%s", it->c_str());
611  FillClusterHist(histname, clust, zv, isPileupEvent, triggers[0]);
612  }
613  }
614  }
615  }
616 
617  PostData(1, fOutput);
618  return true;
619  }
620 
621  //______________________________________________________________________________
622  void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
623  /*
624  * Creating the default pt binning.
625  *
626  * @param binning: Array where to store the results.
627  */
628  std::vector<double> mybinning;
629  std::map<double,double> definitions;
630  definitions.insert(std::pair<double,double>(2.5, 0.1));
631  definitions.insert(std::pair<double,double>(7., 0.25));
632  definitions.insert(std::pair<double,double>(15., 0.5));
633  definitions.insert(std::pair<double,double>(25., 1.));
634  definitions.insert(std::pair<double,double>(40., 2.5));
635  definitions.insert(std::pair<double,double>(50., 5.));
636  definitions.insert(std::pair<double,double>(100., 10.));
637  double currentval = 0;
638  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
639  double limit = id->first, binwidth = id->second;
640  while(currentval < limit){
641  currentval += binwidth;
642  mybinning.push_back(currentval);
643  }
644  }
645  binning.Set(mybinning.size());
646  int ib = 0;
647  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
648  binning[ib++] = *it;
649  }
650 
651  //______________________________________________________________________________
652  void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
653  /*
654  * Creating default z-Vertex binning.
655  *
656  * @param binning: Array where to store the results.
657  */
658  std::vector<double> mybinning;
659  double currentval = -10;
660  mybinning.push_back(currentval);
661  while(currentval <= 10.){
662  currentval += 5.;
663  mybinning.push_back(currentval);
664  }
665  binning.Set(mybinning.size());
666  int ib = 0;
667  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
668  binning[ib++] = *it;
669  }
670 
671  //______________________________________________________________________________
672  void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
673  /*
674  * Creating default z-Vertex binning.
675  *
676  * @param binning: Array where to store the results.
677  */
678  std::vector<double> mybinning;
679  double currentval = -0.8;
680  mybinning.push_back(currentval);
681  while(currentval <= 0.8){
682  currentval += 0.1;
683  mybinning.push_back(currentval);
684  }
685  binning.Set(mybinning.size());
686  int ib = 0;
687  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
688  binning[ib++] = *it;
689  }
690 
691  //______________________________________________________________________________
692  void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
693  const char* title, const TArrayD& binning, const char** labels) {
694  /*
695  * Define an axis with a given binning
696  *
697  * @param axis: Axis to be defined
698  * @param name: Name of the axis
699  * @param title: Title of the axis
700  * @param binning: axis binning
701  * @param labels (@optional): array of bin labels
702  */
703  axis.Set(binning.GetSize()-1, binning.GetArray());
704  axis.SetName(name);
705  axis.SetTitle(title);
706  if(labels){
707  for(int ib = 1; ib <= axis.GetNbins(); ++ib)
708  axis.SetBinLabel(ib, labels[ib-1]);
709  }
710  }
711 
712  //______________________________________________________________________________
713  void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
714  const char* title, int nbins, double min, double max,
715  const char** labels) {
716  /*
717  * Define an axis with number of bins from min to max
718  *
719  * @param axis: Axis to be defined
720  * @param name: Name of the axis
721  * @param title: Title of the axis
722  * @param nbins: Number of bins
723  * @param min: lower limit of the axis
724  * @param max: upper limit of the axis
725  * @param labels (@optional): array of bin labels
726  */
727  axis.Set(nbins, min, max);
728  axis.SetName(name);
729  axis.SetTitle(title);
730  if(labels){
731  for(int ib = 1; ib <= axis.GetNbins(); ++ib)
732  axis.SetBinLabel(ib, labels[ib-1]);
733  }
734  }
735 
736 
737  //______________________________________________________________________________
738  void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
739  double vz, bool isPileup) {
740  /*
741  * Fill event-based histogram
742  *
743  * @param trigger: name of the trigger configuration to be processed
744  * @param vz: z-position of the vertex
745  * @param isPileup: signalises if the event is flagged as pileup event
746  */
747  char histname[1024];
748  sprintf(histname, "hEventHist%s", trigger);
749  try{
750  fHistos->FillTH2(histname, 0., vz);
751  } catch (HistoContainerContentException &e){
752  std::stringstream errormessage;
753  errormessage << "Filling of histogram failed: " << e.what();
754  AliError(errormessage.str().c_str());
755  }
756  if(!isPileup){
757  try{
758  fHistos->FillTH2(histname, 1., vz);
759  } catch(HistoContainerContentException &e){
760  std::stringstream errormessage;
761  errormessage << "Filling of histogram failed: " << e.what();
762  AliError(errormessage.str().c_str());
763  }
764  }
765  }
766 
767  //______________________________________________________________________________
768  void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
769  const AliVTrack* track, double vz, bool isPileup, int cut, bool isMinBias, double jetradius) {
770  /*
771  * Fill track-based histogram with corresponding information
772  *
773  * @param trigger: name of the trigger
774  * @param track: ESD track selected
775  * @param vz: z-position of the vertex
776  * @param isPileup: flag event as pileup event
777  * @param cut: id of the cut (0 = no cut)
778  */
779  double etasign = fSwapEta ? -1. : 1.;
780  double data[7] = {TMath::Abs(track->Pt()), etasign * track->Eta(), track->Phi(), vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
781  double dataMC[7] = {0., 0., 0., vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
782  AliVParticle *assocMC(NULL);
783  if(fMCEvent && (assocMC = fMCEvent->GetTrack(TMath::Abs(track->GetLabel())))){
784  // Select track onl
785  dataMC[0] = TMath::Abs(assocMC->Pt());
786  dataMC[1] = etasign * assocMC->Eta();
787  dataMC[2] = assocMC->Phi();
788  }
789  char histname[1024], histnameAcc[1024], histnameMC[1024], histnameMCAcc[1024];
790  sprintf(histname, "hTrackHist%s", trigger);
791  sprintf(histnameAcc, "hTrackInAcceptanceHist%s", trigger);
792  sprintf(histnameMC, "hMCTrackHist%s", trigger);
793  sprintf(histnameMCAcc, "hMCTrackInAcceptanceHist%s", trigger);
794  if(jetradius > 0.){
795  char *hnames[] = {histname, histnameAcc, histnameMC, histnameMCAcc};
796  for(unsigned int iname = 0; iname < sizeof(hnames)/sizeof(char *); iname++){
797  char *myhname = hnames[iname];
798  sprintf(myhname, "%sRad%02d", myhname, int(jetradius * 10.));
799  }
800  }
801  Bool_t isEMCAL = kFALSE;
802  if(track->IsEMCAL()){
803  // Check if the cluster is matched to only one track
804  AliVCluster *emcclust(NULL);
805  AliDebug(2, Form("cluster id: %d\n", track->GetEMCALcluster()));
806  if(fCaloClusters) {
807  AliDebug(2, "Using calibrated clusters");
808  emcclust = dynamic_cast<AliVCluster *>(fCaloClusters->At(track->GetEMCALcluster()));
809  } else {
810  AliDebug(2, "Using uncalibrated clusters");
811  emcclust = fInputEvent->GetCaloCluster(track->GetEMCALcluster());
812  }
813  if(!emcclust) AliError("Null pointer to EMCal cluster");
814  if(emcclust && emcclust->GetNTracksMatched() <= 1){
815  isEMCAL = kTRUE;
816  }
817  }
818  try{
819  fHistos->FillTHnSparse(histname, data);
820  if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
821  if(isEMCAL){
822  fHistos->FillTHnSparse(histnameAcc, data);
823  if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
824  }
825  } catch (HistoContainerContentException &e){
826  std::stringstream errormessage;
827  errormessage << "Filling of histogram failed: " << e.what();
828  AliError(errormessage.str().c_str());
829  }
830  if(!isPileup){
831  data[4] = 1;
832  dataMC[4] = 1;
833  try{
834  fHistos->FillTHnSparse(histname, data);
835  if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
836  if(isEMCAL){
837  fHistos->FillTHnSparse(histnameAcc, data);
838  if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
839  }
840  } catch (HistoContainerContentException &e){
841  std::stringstream errormessage;
842  errormessage << "Filling of histogram failed: " << e.what();
843  AliError(errormessage.str().c_str());
844  }
845  }
846  }
847 
848  //______________________________________________________________________________
849  void AliAnalysisTaskPtEMCalTrigger::FillClusterHist(const char* histname,
850  const AliVCluster* clust, double vz, bool isPileup, bool isMinBias) {
851  /*
852  * Fill cluster-based histogram with corresponding information
853  *
854  * @param trigger: name of the trigger
855  * @param cluster: the EMCal cluster information
856  * @param vz: z-position of the vertex
857  * @param isPileup: flag event as pileup event
858  */
859  double data[4] = {clust->E(), vz, 0, isMinBias ? 1. : 0.};
860  try{
861  fHistos->FillTHnSparse(histname, data);
862  } catch (HistoContainerContentException &e){
863  std::stringstream errormessage;
864  errormessage << "Filling of histogram failed: " << e.what();
865  AliError(errormessage.str().c_str());
866  }
867  if(!isPileup){
868  data[2] = 1.;
869  try{
870  fHistos->FillTHnSparse(histname, data);
871  } catch (HistoContainerContentException &e){
872  std::stringstream errormessage;
873  errormessage << "Filling of histogram failed: " << e.what();
874  AliError(errormessage.str().c_str());
875  }
876  }
877  }
878 
879  //______________________________________________________________________________
880  void AliAnalysisTaskPtEMCalTrigger::FillMCParticleHist(const char *histname, const AliVParticle * const track, double vz, bool isPileup){
881  /*
882  * Fill histogram for MC-true particles with the information pt, eta and phi
883  *
884  * @param track: the Monte-Carlo track
885  */
886  double data[5] = {TMath::Abs(track->Pt()), track->Eta(), track->Phi(), vz, 0.};
887  fHistos->FillTHnSparse(histname, data);
888  if(!isPileup){
889  data[4] = 1.;
890  fHistos->FillTHnSparse(histname, data);
891  }
892  }
893 
894  //______________________________________________________________________________
895  bool AliAnalysisTaskPtEMCalTrigger::IsTrueTrack(const AliVTrack *const track) const{
896  /*
897  * Check if the track has an associated MC particle, and that the particle is a physical primary
898  * In case of data we do not do the selection at that step (always return true)
899  *
900  * @param track: Track to check
901  * @result: true primary track (true or false)
902  */
903  if(!fMCEvent) return true;
904  AliVParticle *mcassociate = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
905  if(!mcassociate) return false;
906  return fMCEvent->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));
907  }
908 
909  //______________________________________________________________________________
910  void AliAnalysisTaskPtEMCalTrigger::AddESDTrackCuts(AliESDtrackCuts* trackCuts) {
911  /*
912  * Add new track cuts to the task
913  *
914  * @param trackCuts: Object of type AliESDtrackCuts
915  */
916  fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionESD(trackCuts));
917  }
918 
919  //______________________________________________________________________________
920  void AliAnalysisTaskPtEMCalTrigger::AddCutsForAOD(AliESDtrackCuts* trackCuts, UInt_t filterbits) {
921  /*
922  * Add new track cuts to the task
923  *
924  * @param trackCuts: Object of type AliESDtrackCuts
925  */
926  fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionAOD(trackCuts, filterbits));
927  }
928 
929 
930  //______________________________________________________________________________
931  TString AliAnalysisTaskPtEMCalTrigger::BuildTriggerString() {
932  /*
933  * Build trigger string from the trigger maker
934  *
935  * @return: blank-separated string of fired trigger classes
936  */
937  AliDebug(1, "trigger checking");
938  TString result = "";
939  if(HasTriggerType(kJ1)) result += "EJ1 ";
940  if(HasTriggerType(kJ2)) result += "EJ2 ";
941  if(HasTriggerType(kG1)) result += "EG1 ";
942  if(HasTriggerType(kG2)) result += "EG2 ";
943  return result;
944  }
945 
946  //______________________________________________________________________________
947  const AliVVertex* AliAnalysisTaskPtEMCalTrigger::GetSPDVertex() const {
948  /*
949  * Accessor for the SPD vertex, creating transparency for ESDs and AODs
950  *
951  * @return: the spd vertex
952  */
953  AliESDEvent *esd = dynamic_cast<AliESDEvent *>(fInputEvent);
954  if(esd){
955  return esd->GetPrimaryVertexSPD();
956  } else {
957  AliAODEvent *aod = dynamic_cast<AliAODEvent *>(fInputEvent);
958  if(aod){
959  return aod->GetPrimaryVertexSPD();
960  }
961  }
962  return NULL;
963  }
964 
965  //______________________________________________________________________________
966  const AliEmcalJet * AliAnalysisTaskPtEMCalTrigger::FoundTrackInJet(
967  const AliVParticle * const track, AliJetContainer *const jets) const
968  {
969  /*
970  * Correlate track to reconstructed jet
971  *
972  * @param track: particle to be checked
973  * @param jets: container of recontructed jets
974  * @return: The matched jet (NULL if not found)
975  */
976  const AliEmcalJet *result = NULL;
977  const AliEmcalJet *tmp = jets->GetNextAcceptJet(0);
978  while(tmp){
979  if(TrackInJet(track, tmp, jets->GetParticleContainer())){
980  result = tmp;
981  break;
982  }
983  tmp = jets->GetNextAcceptJet();
984  }
985  return result;
986  }
987 
988  //______________________________________________________________________________
989  bool AliAnalysisTaskPtEMCalTrigger::IsInRadius(const AliVParticle *const track, const AliEmcalJet *reconstructedJet, Double_t radius) const {
990  /*
991  * Check if track is in radius around a given jet
992  *
993  * @param track: Track to check
994  * @param reconstructed jet: jet to probe
995  * @param radius: cone radius
996  * @return: result of the test (true if track is inside the cone radius, false otherwise)
997  */
998  return reconstructedJet->DeltaR(track) < radius;
999  }
1000 
1001  //______________________________________________________________________________
1002  bool AliAnalysisTaskPtEMCalTrigger::IsInRadius(const AliVCluster *const clust, const AliEmcalJet *reconstructedJet, Double_t radius) const {
1003  /*
1004  * Check if track is in radius around a given jet
1005  *
1006  * @param track: Track to check
1007  * @param reconstructed jet: jet to probe
1008  * @param radius: cone radius
1009  * @return: result of the test (true if track is inside the cone radius, false otherwise)
1010  */
1011  double vertexPos[3];
1012  fInputEvent->GetPrimaryVertex()->GetXYZ(vertexPos);
1013  TLorentzVector clustVect;
1014  clust->GetMomentum(clustVect, vertexPos);
1015 
1016  Double_t dPhi = reconstructedJet->Phi() - clustVect.Phi();
1017  Double_t dEta = reconstructedJet->Eta() - clustVect.Eta();
1018  dPhi = TVector2::Phi_mpi_pi ( dPhi );
1019 
1020  double deltaR = TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
1021  return deltaR < radius;
1022  }
1023 
1024  //______________________________________________________________________________
1025  bool AliAnalysisTaskPtEMCalTrigger::TrackInJet(const AliVParticle* const track,
1026  const AliEmcalJet* reconstructedJet, const AliParticleContainer * const particles) const {
1027  /*
1028  * Check whether track is among the jet constituents
1029  *
1030  * @param track: track to check
1031  * @param reconstructedJet: reconstructed jet to check
1032  * @param tracks: container with tracks used for jetfinding
1033  * @return: true if found, false otherwise
1034  */
1035  bool found = false;
1036  const AliPicoTrack *picotmp(NULL);
1037  const AliVParticle *tmp(NULL);
1038  for(int ipart = 0; ipart < reconstructedJet->GetNumberOfTracks(); ipart++){
1039  tmp = dynamic_cast<const AliVParticle *>(reconstructedJet->TrackAt(ipart, particles->GetArray()));
1040  if((picotmp = dynamic_cast<const AliPicoTrack *>(tmp))) // handle pico tracks
1041  tmp = picotmp->GetTrack();
1042  if(!tmp->Compare(track)){
1043  found = true;
1044  break;
1045  }
1046  }
1047  return found;
1048  }
1049 
1050  //______________________________________________________________________________
1051  const AliEmcalJet* AliAnalysisTaskPtEMCalTrigger::FoundClusterInJet(const AliVCluster* const clust, AliJetContainer* const jets) const {
1052  /*
1053  * Check whether a cluster is in a radius around a given jet
1054  *
1055  * @param clust: The cluster to check
1056  * @param reconstructedJet: reconstructed jet to check
1057  * @return: the jet containing the cluster (null otherwise)
1058  */
1059  const AliEmcalJet *result = NULL;
1060  const AliEmcalJet *tmp = jets->GetNextAcceptJet(0);
1061  while(tmp){
1062  if(ClusterInJet(clust, tmp, jets->GetClusterContainer())){
1063  result = tmp;
1064  break;
1065  }
1066  tmp =jets->GetNextAcceptJet();
1067  }
1068  return result;
1069  }
1070 
1071  //______________________________________________________________________________
1072  bool AliAnalysisTaskPtEMCalTrigger::ClusterInJet(const AliVCluster* const clust,
1073  const AliEmcalJet* reconstructedJet, const AliClusterContainer* const clusters) const {
1074  /*
1075  * Check whether cluster is among the jet constituents
1076  *
1077  * @param track: track to check
1078  * @param reconstructedJet: reconstructed jet to check
1079  * @param clusters: the cluster container
1080  * @return: true if found, false otherwise
1081  */
1082  bool found = false;
1083  const AliVCluster *tmp(NULL);
1084  for(int ipart = 0; ipart < reconstructedJet->GetNumberOfTracks(); ipart++){
1085  tmp = dynamic_cast<const AliVCluster *>(reconstructedJet->ClusterAt(ipart, clusters->GetArray()));
1086  if(!tmp->Compare(clust)){
1087  found = true;
1088  break;
1089  }
1090  }
1091  return found;
1092  }
1093 
1094  //______________________________________________________________________________
1095  void EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger::AddJetContainerName(const Char_t* contname, Bool_t isMC) {
1096  /*
1097  * Add new Jet input container to the analysis task
1098  *
1099  * @param contname: Name of the container
1100  * @param isMC: Defines whether the container is for MC truth or not
1101  */
1102  TList &mycontainer = isMC ? fJetContainersMC : fJetContainersData;
1103  mycontainer.Add(new TObjString(contname));
1104  }
1105 }
title
const char * title
Definition: MakeQAPdf.C:26
AliEMCalPtTaskTrackSelectionAOD.h
Implement virtual track selection for AOD analysis.
AliAnalysisTaskPtEMCalTrigger.h
EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger::AddJetContainerName
void AddJetContainerName(const Char_t *contname, Bool_t isMC=kFALSE)
AliEMCalHistoContainer.h
Declarartion of class AliEMCalHistoContainer.
AliEMCalPtTaskTrackSelectionESD.h
Declaration of class AliEMCalPtTaskTrackSelectionESD.
AliEMCalPtTaskVTrackSelection.h
Declartion of class AliEMCalPtTaskVTrackSelection.
ClassImp
ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger) namespace EMCalTriggerPtAnalysis
Definition: AliAnalysisTaskPtEMCalTrigger.cxx:65