AliPhysics  8d00e07 (8d00e07)
AddTaskCentralTracks.C
Go to the documentation of this file.
1 
12 #include <AliAnalysisTaskSE.h>
13 #include <AliESDtrackCuts.h>
14 class TList;
15 class TH2D;
16 class TH1D;
17 
85 {
86 public:
87  enum {
88  kValidTrigger = 0x1,
90  };
103  CentralMultTask(const char* name,
104  Double_t maxEta=2,
105  Double_t maxVtx=10);
110  virtual ~CentralMultTask();
116  void SetUseTracklets(Bool_t use) { fUseTracklets = use; }
121  virtual void Init() {}
127  virtual void UserCreateOutputObjects();
133  virtual void UserExec(Option_t* option);
141  virtual void Terminate(Option_t* option);
142 protected:
151  UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
152 
153  TH2D* fHist; // Per-event d2n/deta/dphi
154  TH1D* fAll; // Summed d2n/deta/dphi - all track(let)s
155  TH1D* fGlobal; // Summed d2n/deta/dphi - global tracks
156  TH1D* fITS; // Summed d2n/deta/dphi - ITS tracks
157  TH1D* fTracklets; // Summed d2n/deta/dphi - SPD tracklets
158  TH1D* fNEventsTr; // Number of triggered events per vertex bin
159  TH1D* fNEventsVtx;// Number of triggered+vertex events per vertex
160  TList* fOutput; // Output list
161  AliESDtrackCuts fQGlo; // Quality cut on ITS+TPC
162  AliESDtrackCuts fQITS; // Quality cut on ITS standalone (not ITS+TPC)
163  AliESDtrackCuts fDCAwSPD; // DCA for traks with SPD hits
164  AliESDtrackCuts fDCAwoSPD; // DCA for traks without SPD hits
165  AliESDtrackCuts fIsPrimary; // Primary particle
166  Bool_t fUseTracklets; // Whether to count tracklets or not
167  Bool_t fDebug; // Debug flag
168 
169  ClassDef(CentralMultTask,1); // Determine multiplicity in central area
170 };
171 
172 //====================================================================
173 #include <TMath.h>
174 #include <TH2D.h>
175 #include <TH1D.h>
176 #include <THStack.h>
177 #include <TList.h>
178 #include <AliAnalysisManager.h>
179 #include <AliESDEvent.h>
180 #include <AliAODHandler.h>
181 #include <AliAODEvent.h>
182 #include <AliESDInputHandler.h>
183 #include <AliESDtrack.h>
184 #include <AliMultiplicity.h>
185 
186 //____________________________________________________________________
188  : AliAnalysisTaskSE(),
189  fHist(0),
190  fAll(0),
191  fGlobal(0),
192  fITS(0),
193  fTracklets(0),
194  fNEventsTr(0),
195  fNEventsVtx(0),
196  fOutput(0),
197  fUseTracklets(false),
198  fDebug(false)
199 {}
200 
201 //____________________________________________________________________
202 inline CentralMultTask::CentralMultTask(const char* /* name */,
203  Double_t maxEta,
204  Double_t maxVtx)
205  : AliAnalysisTaskSE("Central"),
206  fHist(0),
207  fAll(0),
208  fGlobal(0),
209  fITS(0),
210  fTracklets(0),
211  fNEventsTr(0),
212  fNEventsVtx(0),
213  fOutput(0),
214  fUseTracklets(true),
215  fDebug(false)
216 {
217  Int_t nBins = 40;
218  fHist = new TH2D("Central", "d^{2}N_{ch}/d#etad#phi in central region",
219  nBins, -maxEta, maxEta, 20, 0, 2*TMath::Pi());
220  fHist->SetDirectory(0);
221  fHist->SetXTitle("#eta");
222  fHist->SetYTitle("#phi [radians]");
223  fHist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
224  fHist->SetStats(0);
225  fHist->Sumw2();
226 
227  fAll = new TH1D("all", "Central region", nBins, -maxEta, maxEta);
228  fAll->SetDirectory(0);
229  fAll->SetXTitle("#eta");
230  fAll->SetYTitle("dN_{ch}/d#eta");
231  fAll->Sumw2();
232  fAll->SetFillColor(kGray);
233  fAll->SetFillStyle(3001);
234  fAll->SetMarkerStyle(28);
235  fAll->SetMarkerColor(kGray);
236  fAll->SetStats(0);
237 
238  fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
239  fGlobal->SetTitle("Global tracks");
240  fGlobal->SetDirectory(0);
241  fGlobal->Sumw2();
242  fGlobal->SetFillColor(kRed+1);
243  fGlobal->SetFillStyle(3001);
244  fGlobal->SetMarkerStyle(28);
245  fGlobal->SetMarkerColor(kRed+1);
246  fGlobal->SetStats(0);
247 
248  fITS = static_cast<TH1D*>(fAll->Clone("its"));
249  fITS->SetTitle("ITS tracks");
250  fITS->SetDirectory(0);
251  fITS->Sumw2();
252  fITS->SetFillColor(kGreen+1);
253  fITS->SetFillStyle(3001);
254  fITS->SetMarkerStyle(28);
255  fITS->SetMarkerColor(kGreen+1);
256  fITS->SetStats(0);
257 
258  fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
259  fTracklets->SetTitle("SPD tracklets");
260  fTracklets->SetDirectory(0);
261  fTracklets->Sumw2();
262  fTracklets->SetFillColor(kBlue+1);
263  fTracklets->SetFillStyle(3001);
264  fTracklets->SetMarkerStyle(28);
265  fTracklets->SetMarkerColor(kBlue+1);
266  fTracklets->SetStats(0);
267 
268  fNEventsTr = new TH1D("nEventsTr",
269  "Events per vertex bin", 10, -maxVtx, maxVtx);
270  fNEventsTr->SetXTitle("v_{z}");
271  fNEventsTr->SetYTitle("Events");
272  fNEventsTr->SetDirectory(0);
273 
274  fNEventsVtx = new TH1D("nEventsVtx",
275  "Events per vertex bin", 10, -maxVtx, maxVtx);
276  fNEventsVtx->SetXTitle("v_{z}");
277  fNEventsVtx->SetYTitle("Events");
278  fNEventsVtx->SetDirectory(0);
279 
280  // --- Global (ITS+TPC) track quality cuts
281  // TPC
282  fQGlo.SetMinNClustersTPC(70);
283  fQGlo.SetMaxChi2PerClusterTPC(4);
284  fQGlo.SetAcceptKinkDaughters(kFALSE);
285  fQGlo.SetRequireTPCRefit(kTRUE);
286  // ITS
287  fQGlo.SetRequireITSRefit(kTRUE);
288  fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
289  fQGlo.SetEtaRange(-maxEta, maxEta);
290 
291  // --- ITS standalone track quality cuts
292  fQITS.SetRequireITSRefit(kTRUE);
293  fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
294  fQITS.SetEtaRange(-maxEta, maxEta);
295 
296  // -- Distance-of-Closest-Approach cuts for tracks w/ITS hits
297  fDCAwSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
298  AliESDtrackCuts::kAny);
299  fDCAwSPD.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
300  fDCAwSPD.SetMaxDCAToVertexZ(0.5);
301  fDCAwSPD.SetEtaRange(-maxEta, maxEta);
302 
303  // -- Distance-of-Closest-Approach cuts for tracks w/o ITS hits
304  fDCAwoSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
305  AliESDtrackCuts::kNone);
306  fDCAwoSPD.SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
307  fDCAwoSPD.SetMaxDCAToVertexZ(0.5);
308  fDCAwoSPD.SetEtaRange(-maxEta, maxEta);
309 
310  // -- Primary track cut
311  // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
312  // Quality
313  fIsPrimary.SetMinNClustersTPC(70);
314  fIsPrimary.SetMaxChi2PerClusterTPC(4);
315  fIsPrimary.SetAcceptKinkDaughters(kFALSE);
316  fIsPrimary.SetRequireTPCRefit(kTRUE);
317  fIsPrimary.SetRequireITSRefit(kTRUE);
318  fIsPrimary.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
319  AliESDtrackCuts::kAny);
320  // Dca:
321  fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
322  fIsPrimary.SetMaxDCAToVertexZ(2);
323  fIsPrimary.SetDCAToVertex2D(kFALSE);
324  fIsPrimary.SetRequireSigmaToVertex(kFALSE);
325 
326  // Output slot #1 writes into a TH1 container
327  DefineOutput(1, TList::Class());
328 }
329 
330 //____________________________________________________________________
332 {
333  if (fHist) delete fHist;
334  if (fAll) delete fAll;
335  if (fGlobal) delete fGlobal;
336  if (fITS) delete fITS;
337  if (fTracklets) delete fTracklets;
338 }
339 
340 //________________________________________________________________________
341 inline void
343 {
344  // Create histograms
345  // Called once (on the worker node)
346 
347  fOutput = new TList;
348  fOutput->SetName(GetName());
349  fOutput->SetOwner();
350 
351  fOutput->Add(fAll);
352  fOutput->Add(fGlobal);
353  fOutput->Add(fITS);
354  fOutput->Add(fTracklets);
355  fOutput->Add(fNEventsTr);
356  fOutput->Add(fNEventsVtx);
357 
358  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
359  AliAODHandler* ah =
360  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
361  if (!ah) AliFatal("No AOD output handler set in analysis manager");
362 
363  ah->AddBranch("TH2D", &fHist);
364 
365  // Post data for ALL output slots >0 here, to get at least an empty histogram
366  PostData(1, fOutput);
367 }
368 
369 //____________________________________________________________________
370 inline UShort_t
372 {
373  // Do some fast cuts first
374  vz = 0;
375 
376  // Get the analysis manager - should always be there
377  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
378  if (!am) {
379  AliWarning("No analysis manager defined!");
380  return 0;
381  }
382 
383  // Get the input handler - should always be there
384  AliInputEventHandler* ih =
385  static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
386  if (!ih) {
387  AliWarning("No input handler");
388  return 0;
389  }
390 
391  // Trigger mask
392  UInt_t mask = ih->IsEventSelected();
393  Bool_t isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
394  UShort_t ret = 0;
395  if (isMinBias) ret |= kValidTrigger;
396 
397  // check for good reconstructed vertex
398  if (!(esd.GetPrimaryVertex()->GetStatus())) {
399  if (fDebug)
400  AliWarning("Primary vertex has bad status");
401  return ret;
402  }
403  if (!(esd.GetPrimaryVertexSPD()->GetStatus())) {
404  if (fDebug)
405  AliWarning("Primary SPD vertex has bad status");
406  return ret;
407  }
408 
409  // if vertex is from spd vertexZ, require more stringent cut
410  if (esd.GetPrimaryVertex()->IsFromVertexerZ()) {
411  if (esd.GetPrimaryVertex()->GetDispersion() > 0.02 ||
412  esd.GetPrimaryVertex()->GetZRes() > 0.25) {
413  if (fDebug)
414  AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)",
415  esd.GetPrimaryVertex()->GetDispersion(),
416  esd.GetPrimaryVertex()->GetZRes()));
417  return ret;
418  }
419  }
420  // One can add a cut on the vertex z position here
421  vz = esd.GetPrimaryVertex()->GetZ();
422  Double_t vl = fNEventsVtx->GetXaxis()->GetXmin();
423  Double_t vh = fNEventsVtx->GetXaxis()->GetXmax();
424  if (vz < vl || vz > vh) {
425  if (fDebug)
426  AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
427  return ret;
428  }
429  ret |= kValidVertex;
430 
431  return ret;
432 }
433 
434 //____________________________________________________________________
435 inline void
437 {
438  // Main loop
439  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
440  if (!esd) {
441  AliError("Cannot get the ESD event");
442  return;
443  }
444 
445  fHist->Reset();
446 
447  // Check event
448  Double_t vz = 0;
449  UShort_t eventFlags = CheckEvent(*esd, vz);
450  if (!(eventFlags & kValidTrigger))
451  return; // No trigger
452  else {
453  fNEventsTr->Fill(vz);
454  if (eventFlags & kValidVertex)
455  fNEventsVtx->Fill(vz);
456  else
457  return; // No trigger or no vertex
458  }
459 
460 
461  // flags for secondary and rejected tracks
462  // set this bit in ESD tracks if it is rejected by a cut
463  const int kRejBit = BIT(15);
464  // set this bit in ESD tracks if it is secondary according to a cut
465  const int kSecBit = BIT(16);
466 
467  Int_t total = 0;
468  Int_t nESDTracks = esd->GetNumberOfTracks();
469  // Loop over ESD tracks
470  for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
471 
472  AliESDtrack* track = esd->GetTrack(i);
473 
474  // if track is a secondary from a V0, flag as a secondary
475  if (track->IsOn(AliESDtrack::kMultInV0)) {
476  track->SetBit(kSecBit);
477  continue;
478  }
479 
480  Double_t eta = track->Eta();
481  Double_t phi = track->Phi();
482  // check tracks with ITS part
483  if (track->IsOn(AliESDtrack::kITSin)) {
484 
485  // Check ITS pure stand-alone - and if it is, reject it
486  if (track->IsOn(AliESDtrack::kITSpureSA)) {
487  track->SetBit(kRejBit);
488  continue;
489  }
490 
491  // Check DCA - if not close enough reject it as secondary
492  if (!fDCAwSPD.AcceptTrack(track) &&
493  !fDCAwoSPD.AcceptTrack(track)) {
494  track->SetBit(kSecBit);
495  continue;
496  }
497 
498  // Check if this is an ITS complementary - no TPC in
499  if (!track->IsOn(AliESDtrack::kTPCin)) {
500  if (fQITS.AcceptTrack(track)) { // Check ITS quality
501  fITS->Fill(eta);
502  fAll->Fill(eta);
503  }
504  else
505  track->SetBit(kRejBit);
506  }
507  else { // Not ITS SA or TPC+ITS
508  if (fQGlo.AcceptTrack(track)) { // Check ITS quality
509  fGlobal->Fill(eta);
510  fAll->Fill(eta);
511  }
512  else
513  track->SetBit(kRejBit);
514  }
515  if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
516 
517  // fill d2n/detadphi
518  fHist->Fill(eta, phi);
519  }
520  }
521 
522  // Get multiplicity from SPD tracklets
523  const AliMultiplicity* spdmult = esd->GetMultiplicity();
524  for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i){
525  // if counting tracks+tracklets,
526  // check if clusters were already used in tracks
527  Int_t id1,id2;
528  spdmult->GetTrackletTrackIDs(i,0,id1,id2);
529  AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
530  AliESDtrack* tr2 = id2>=0 ? esd->GetTrack(id2) : 0;
531  //
532  if ((tr1 && tr1->TestBit(kSecBit)) || // Flagged as secondary
533  (tr2 && tr2->TestBit(kSecBit)) || // Flagged as secondary
534  (tr1 && !tr1->TestBit(kRejBit)) || // already accounted for
535  (tr2 && !tr2->TestBit(kRejBit))) // already accounted for
536  continue;
537  ++total;
538  Double_t eta = spdmult->GetEta(i);
539  Double_t phi = spdmult->GetPhi(i);
540 
541  // Increment d2n/detadphi from an SPD tracklet
542  if (fUseTracklets) {
543  fHist->Fill(eta, phi);
544  fAll->Fill(eta);
545  total++;
546  }
547  fTracklets->Fill(eta);
548  }
549  if (fDebug) AliInfo(Form("A total of %d tracks", total));
550 
551  // NEW HISTO should be filled before this point, as PostData puts the
552  // information for this iteration of the UserExec in the container
553  PostData(1, fOutput);
554 }
555 
556 
557 //________________________________________________________________________
558 inline void
560 {
561  // Draw result to screen, or perform fitting, normalizations Called
562  // once at the end of the query
563 
564  fOutput = dynamic_cast<TList*> (GetOutputData(1));
565  if(!fOutput) {
566  AliError("Could not retrieve TList fOutput");
567  return;
568  }
569 
570  TH1D* all = static_cast<TH1D*>(fOutput->FindObject("all"));
571  TH1D* global = static_cast<TH1D*>(fOutput->FindObject("global"));
572  TH1D* its = static_cast<TH1D*>(fOutput->FindObject("its"));
573  TH1D* tracklets = static_cast<TH1D*>(fOutput->FindObject("tracklets"));
574  TH1D* eventsTr = static_cast<TH1D*>(fOutput->FindObject("nEventsTr"));
575  TH1D* eventsVtx = static_cast<TH1D*>(fOutput->FindObject("nEventsVtx"));
576 
577  Int_t nTriggers = eventsTr->GetEntries();
578  Int_t nVertex = eventsVtx->GetEntries();
579  if (nTriggers <= 0 || nVertex <= 0) {
580  AliWarning("No data in the events histogram!");
581  nTriggers = 1;
582  }
583 
584  all ->Scale(1. / nTriggers, "width");
585  global ->Scale(1. / nTriggers, "width");
586  its ->Scale(1. / nTriggers, "width");
587  tracklets->Scale(1. / nTriggers, "width");
588 
589  THStack* stack = new THStack("components", "Components");
590  if (fUseTracklets) stack->Add(tracklets);
591  stack->Add(global);
592  stack->Add(its);
593 
594  fOutput->Add(stack);
595 }
596 
597 //========================================================================
604 inline AliAnalysisTask*
606 {
607  // analysis manager
608  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
609 
610  // Make our object. 2nd argumenent is absolute max Eta
611  // 3rd argument is absolute max Vz
612  CentralMultTask* task = new CentralMultTask("Global", 2, 10);
613  // if physics selection performed in UserExec(),
614  // this line should be commented
615  // task->SelectCollisionCandidates(AliVEvent::kMB);
616  mgr->AddTask(task);
617 
618  // create containers for input/output
619  AliAnalysisDataContainer *output =
620  mgr->CreateContainer("Central", TList::Class(),
621  AliAnalysisManager::kOutputContainer,
622  AliAnalysisManager::GetCommonFileName());
623 
624  // connect input/output
625  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
626  mgr->ConnectOutput(task, 1, output);
627 
628  return task;
629 }
630 
631 
632 //________________________________________________________________________
633 //
634 // EOF
635 //
double Double_t
Definition: External.C:58
virtual void Terminate(Option_t *option)
AliESDtrackCuts fDCAwSPD
AliStack * stack
virtual void UserExec(Option_t *option)
AliESDtrackCuts fDCAwoSPD
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
AliESDtrackCuts fIsPrimary
UShort_t CheckEvent(const AliESDEvent &esd, Double_t &vz)
void SetUseTracklets(Bool_t use)
Definition: External.C:228
Definition: External.C:212
virtual void UserCreateOutputObjects()
AliAnalysisTask * AddTaskCentralTracks()
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
AliESDtrackCuts fQITS
bool Bool_t
Definition: External.C:53
AliESDtrackCuts fQGlo