AliPhysics  master (3d17d9d)
AliEmcalTRDTrackingTask.cxx
Go to the documentation of this file.
1 #include <cstring>
2 #include <iostream>
3 
4 #include <THnSparse.h>
5 #include <TMath.h>
6 #include <TString.h>
7 
8 #include <AliESDtrack.h>
9 #include <AliAODTrack.h>
10 #include <AliAODMCParticle.h>
11 #include <AliLog.h>
12 #include <AliAnalysisManager.h>
13 #include <AliVEventHandler.h>
14 
15 #include "AliMCParticleContainer.h"
16 
18 
20 
21 
25  AliAnalysisTaskEmcalLight("AliEmcalTrackingQA", kTRUE),
26  fIsEsd(kFALSE),
27  fGeneratorLevel(nullptr),
28  fDetectorLevel(nullptr),
29  fPtHistBins(),
30  fPtResHistBins(),
31  fIntegerHistBins(),
32  fChargeHistBins(),
33  fTracks(nullptr)
34 {
35  SetMakeGeneralHistograms(kTRUE);
36 
37  GenerateHistoBins();
38 }
39 
44  AliAnalysisTaskEmcalLight("AliEmcalTrackingQA", kTRUE),
45  fIsEsd(kFALSE),
46  fGeneratorLevel(nullptr),
47  fDetectorLevel(nullptr),
48  fPtHistBins(),
49  fPtResHistBins(),
50  fIntegerHistBins(),
51  fChargeHistBins(),
52  fTracks(nullptr)
53 {
55 
57 }
58 
63 {
64 }
65 
70 {
71  GenerateFixedBinArray(6, 0.0, 0.3, fPtHistBins, false);
72  GenerateFixedBinArray(7, 0.3, 1.0, fPtHistBins, false);
73  GenerateFixedBinArray(10, 1.0, 3.0, fPtHistBins, false);
74  GenerateFixedBinArray(14, 3.0, 10.0, fPtHistBins, false);
75  GenerateFixedBinArray(10, 10.0, 20.0, fPtHistBins, false);
76  GenerateFixedBinArray(15, 20.0, 50.0, fPtHistBins, false);
77  GenerateFixedBinArray(40, 50.0, 250.0, fPtHistBins, false);
78  GenerateFixedBinArray(10, 250.0, 350.0, fPtHistBins);
79 
80 
81  GenerateFixedBinArray(50, 0.00, 0.05, fPtResHistBins, false);
82  GenerateFixedBinArray(25, 0.05, 0.10, fPtResHistBins, false);
83  GenerateFixedBinArray(25, 0.10, 0.20, fPtResHistBins, false);
84  GenerateFixedBinArray(30, 0.20, 0.50, fPtResHistBins, false);
85  GenerateFixedBinArray(25, 0.50, 1.00, fPtResHistBins, false);
86  GenerateFixedBinArray(20, 1.00, 2.00, fPtResHistBins);
87 
89 
90  GenerateFixedBinArray(2, -1.1, 1.1, fChargeHistBins, true);
91 }
92 
97 {
99 
100  if (fParticleCollArray.empty()) {
101  AliErrorStream() << "This task needs at least one particle container! The task won't run." << std::endl;
102  fInhibit = kTRUE;
103  return;
104  }
105 
106  AliDebugStream(3) << "Loading the detector track container" << std::endl;
107  if (!fDetectorLevel) {
108  auto iter = fParticleCollArray.find("detector");
109  if (iter == fParticleCollArray.end()) {
110  AliErrorStream() << "This task needs at least one particle container named 'detector'! The task won't run." << std::endl;
111  fInhibit = kTRUE;
112  return;
113  }
114  else {
115  fDetectorLevel = static_cast<AliTrackContainer*>(iter->second);
116  }
117  }
118 
119  AliDebugStream(3) << "Loading the generator particle container" << std::endl;
120  if (!fGeneratorLevel) {
121  auto iter = fParticleCollArray.find("generator");
122  if (iter == fParticleCollArray.end()) {
123  AliInfoStream() << "No particle container named 'generator' was found. Assuming this is not a MC production." << std::endl;
124  }
125  else {
126  fGeneratorLevel = static_cast<AliMCParticleContainer*>(iter->second);
127  }
128  }
129 
130  AliDebugStream(3) << "Allocating histograms" << std::endl;
132 }
133 
135 {
137  if (!fDetectorLevel->GetArray()) {
138  AliErrorStream() << "Could not load track array! The task won't run." << std::endl;
139  fInhibit = kTRUE;
140  return;
141  }
142  if (fDetectorLevel->GetArray()->GetClass()->InheritsFrom("AliESDtrack")) fIsEsd = kTRUE;
143 }
144 
151 THnSparse* AliEmcalTRDTrackingTask::GenerateTHnSparse(const char* name, const std::vector<std::tuple<std::string, std::vector<Double_t>::iterator, std::vector<Double_t>::iterator>>& axis)
152 {
153  std::vector<int> nbins;
154  for (auto a : axis) nbins.push_back(int(std::get<2>(a) - std::get<1>(a) - 1));
155 
156  THnSparse* h = new THnSparseF(name, name, nbins.size(), &nbins[0]);
157  Int_t i = 0;
158  for (auto a : axis) {
159  h->GetAxis(i)->SetTitle(std::get<0>(a).c_str());
160  h->SetBinEdges(i, &(*(std::get<1>(a))));
161  i++;
162  }
163 
164  return h;
165 }
166 
171 {
172  typedef std::vector<Double_t>::iterator my_iterator;
173 
174  std::vector<std::tuple<std::string, my_iterator, my_iterator>> axis;
175 
176 
177  axis.push_back(std::make_tuple("#it{p}_{T} (GeV/#it{c})", fPtHistBins.begin(), fPtHistBins.end()));
178  axis.push_back(std::make_tuple("track type", fIntegerHistBins.begin(), fIntegerHistBins.begin() + 5));
179  axis.push_back(std::make_tuple("number of tracklets", fIntegerHistBins.begin(), fIntegerHistBins.begin()+8));
180  axis.push_back(std::make_tuple("#sigma(#it{p}_{T}) / #it{p}_{T}", fPtResHistBins.begin(), fPtResHistBins.end()));
181  axis.push_back(std::make_tuple("charge", fChargeHistBins.begin(), fChargeHistBins.end()));
182 
183  fTracks = GenerateTHnSparse("fTracks", axis);
184 
185  fOutput->Add(fTracks);
186 }
187 
191 void AliEmcalTRDTrackingTask::FillDetectorLevelTHnSparse(Double_t trackPt, Double_t sigma1OverPt, Byte_t trackType, Int_t ntracklets, Double_t trackCharge)
192 {
193  AliDebugStream(10) << "Filling detector level THnSparse" << std::endl;
194  std::vector<Double_t> contents(fTracks->GetNdimensions());
195 
196  for (Int_t i = 0; i < fTracks->GetNdimensions(); i++) {
197  TString title(fTracks->GetAxis(i)->GetTitle());
198  if (title=="#it{p}_{T} (GeV/#it{c})")
199  contents[i] = trackPt;
200  else if (title=="#sigma(#it{p}_{T}) / #it{p}_{T}")
201  contents[i] = sigma1OverPt*trackPt;
202  else if (title=="track type")
203  contents[i] = trackType;
204  else if (title == "number of tracklets")
205  contents[i] = ntracklets;
206  else if (title=="charge")
207  contents[i] = trackCharge;
208  else
209  AliWarning(Form("Unable to fill dimension %s of histogram %s!", title.Data(), fTracks->GetName()));
210  }
211 
212  fTracks->Fill(&contents[0]);
213 }
214 
219 {
220  AliDebugStream(1) << "Called: Tracks [" << fDetectorLevel->GetNTracks() << "], Accepted [" << fDetectorLevel->GetNAcceptedTracks() << "]\n";
221  auto iterable = fDetectorLevel->accepted_momentum();
222  for (auto trackIterator = iterable.begin(); trackIterator != iterable.end(); trackIterator++) {
223  auto track = trackIterator->second;
224  Byte_t type = fDetectorLevel->GetTrackType(track);
225  AliDebugStream(2) << "Next track of type " << static_cast<Int_t>(type) << std::endl;
226  Byte_t ntracklets = 0;
227  if (type <= 3) {
228  Double_t sigma = 0;
229 
230  if (fIsEsd) {
231  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
232  if (esdTrack){
233  sigma = TMath::Sqrt(esdTrack->GetSigma1Pt2());
234  ntracklets = esdTrack->GetTRDntracklets();
235  }
236  }
237  else { // AOD
238  AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(track);
239  if(!aodtrack) AliFatal("Not a standard AOD");
240 
241  AliExternalTrackParam exParam;
242 
243  //get covariance matrix
244  Double_t cov[21] = {0,};
245  aodtrack->GetCovMatrix(cov);
246  Double_t pxpypz[3] = {0,};
247  aodtrack->PxPyPz(pxpypz);
248  Double_t xyz[3] = {0,};
249  aodtrack->GetXYZ(xyz);
250  Short_t sign = aodtrack->Charge();
251  exParam.Set(xyz,pxpypz,cov,sign);
252 
253  sigma = TMath::Sqrt(exParam.GetSigma1Pt2());
254  ntracklets = track->GetTRDntrackletsPID();
255  }
256 
257  // Gold condition:
258  // - at least 3 TRD tracklets (with this cut track without TRD in global track fit is at % level)
259  if(!(track->GetStatus() & AliVTrack::kTRDupdate)) type += 2;
260 
261  Int_t label = TMath::Abs(track->GetLabel());
262  AliDebugStream(10) << "Track " << trackIterator.current_index() << " with type " << int(type) << " and label " << label <<
263  ", pt = " << track->Pt() << std::endl;
264 
265  double pt = track->Pt();
266 
267  FillDetectorLevelTHnSparse(pt, sigma, type, ntracklets, track->Charge());
268  }
269  else {
270  AliErrorStream() << "Track " << trackIterator.current_index() << " has type " << type << " not recognized!" << std::endl;
271  }
272  }
273 
274  return kTRUE;
275 }
276 
283 {
284  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
285  if (!mgr) {
286  AliErrorClassStream() << "No analysis manager to connect to." << std::endl;
287  return nullptr;
288  }
289 
290  // Check the analysis type using the event handlers connected to the analysis manager
291  AliVEventHandler* handler = mgr->GetInputEventHandler();
292  if (!handler) {
293  AliErrorClassStream() << "This task requires an input event handler" << std::endl;
294  return nullptr;
295  }
296  else {
297  AliInfoClassStream() << "Event handler ok!" << std::endl;
298  }
299 
301 
302  if (handler->InheritsFrom("AliESDInputHandler")) {
303  dataType = kESD;
304  AliInfoClassStream() << "Data type is ESD." << std::endl;
305  }
306  else if (handler->InheritsFrom("AliAODInputHandler")) {
307  dataType = kAOD;
308  AliInfoClassStream() << "Data type is AOD." << std::endl;
309  }
310 
311  TString track_name("tracks");
312  if (dataType == kESD) track_name = "Tracks";
313 
314  // Init the task and do settings
315  TString name("AliEmcalTRDTrackingTask");
316  if(strlen(suffix)) name += TString::Format("_%s", suffix);
317  AliInfoClassStream() << "Allocating task." << std::endl;
319  qaTask->SetVzRange(-10,10);
320  AliInfoClassStream() << "Task allocated, setting containers." << std::endl;
321  qaTask->AddParticleContainer(track_name.Data(), "detector");
322 
323  AliInfoClassStream() << "Containers ok, adding task to the analysis manager." << std::endl;
324  // Final settings, pass to manager and set the containers
325  mgr->AddTask(qaTask);
326 
327  AliInfoClassStream() << "Task added, setting input/output." << std::endl;
328  // Create containers for input/output
329  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
330 
331  TString contName(name);
332  contName += "_histos";
333  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(),
334  TList::Class(),AliAnalysisManager::kOutputContainer,
335  Form("%s", AliAnalysisManager::GetCommonFileName()));
336  mgr->ConnectInput (qaTask, 0, cinput1 );
337  mgr->ConnectOutput (qaTask, 1, coutput1 );
338 
339 
340  AliInfoClassStream() << "Task configuration done." << std::endl;
341  return qaTask;
342 }
343 
std::vector< Double_t > fChargeHistBins
Bool_t fInhibit
!inhibit execution of the task
double Double_t
Definition: External.C:58
EDataType_t
Switch for the data type.
Int_t GetNTracks() const
const char * title
Definition: MakeQAPdf.C:27
Container with name, TClonesArray and cuts for particles.
std::vector< Double_t > fPtHistBins
! pt bins
THnSparse * fTracks
! charge bins
static std::vector< double > GenerateFixedBinArray(int n, double min, double max, bool last=true)
AliParticleContainer * AddParticleContainer(EMCAL_STRINGVIEW branchName, EMCAL_STRINGVIEW contName="")
Double_t * sigma
std::vector< Double_t > fIntegerHistBins
! integer bins
int Int_t
Definition: External.C:63
void FillDetectorLevelTHnSparse(Double_t trackPt, Double_t sigma1OverPt, Byte_t trackType, Int_t ntracklets, Double_t trackCharge)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
TClonesArray * GetArray() const
THnSparse * GenerateTHnSparse(const char *name, const std::vector< std::tuple< std::string, std::vector< Double_t >::iterator, std::vector< Double_t >::iterator >> &axis)
std::vector< Double_t > fPtResHistBins
! pt res bins
AliTrackContainer * fDetectorLevel
! detector level container
short Short_t
Definition: External.C:23
Char_t GetTrackType(const AliVTrack *track) const
static AliEmcalTRDTrackingTask * AddTaskTRDTracking(const char *suffix="")
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
void SetVzRange(Double_t min, Double_t max)
const AliTrackIterableMomentumContainer accepted_momentum() const
const Int_t nbins
bool Bool_t
Definition: External.C:53
AliMCParticleContainer * fGeneratorLevel
! generator level container
Container for MC-true particles within the EMCAL framework.
Bool_t fIsEsd
! whether it is ESD data