2 * @file AddTaskCentralTracks.C
3 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
4 * @date Fri Jan 28 10:22:26 2011
6 * @brief Class and script to add a multiplicity task for the central
9 * @ingroup pwg2_forward_scripts_tasks
12 #include <AliAnalysisTaskSE.h>
13 #include <AliESDtrackCuts.h>
19 * Task to determine the
21 * \left.\frac{d^2N_{ch}}{d\eta d\phi}\right|_{central}
23 * from tracks and tracklets.
25 * First, global tracks are investigated. The requirements on global
27 * - @f$ n_{TPC clusters} \ge 70@f$
28 * - @f$ \chi^2_{TPC cluster} \le 4@f$
30 * - Re-fit of TPC tracks
31 * - Re-fit of ITS tracks
32 * - No requirement on SPD clusters
34 * Secondly, ITS stand-alone tracks are investigated. The
35 * requirements on ITS standalone tracks are
36 * - Re-fit of ITS tracks
37 * - No requirement on SPD clusters
39 * Tracks that does not meet these quality conditions are flagged as
42 * Both kinds of tracks (global and ITS standalone) have requirements
43 * on the distance of closest approach (DCA) to the interaction
44 * vertex @f$(v_x,v_y,v_z)@f$.
46 * For tracks with SPD clusters:
47 * - @f$ DCA_{xy} < 0.0182+0.0350/p_t^{1.01}@f$
48 * - @f$ DCA_{z} < 0.5@f$
50 * For tracks without SPD clusters
51 * - @f$ DCA_{xy} < 1.5(0.0182+0.0350/p_t^{1.01})@f$
52 * - @f$ DCA_{z} < 0.5@f$
54 * Tracks that does not meet these DCA requirements are flagged as
57 * Thirdly, the number of SPD tracklets are investigated. If the
58 * tracklet is associated with a track, and that track has already
59 * been used, then that tracklet is ignored
61 * An @f$(\eta,\phi)@f$ per-event histogram is then filled from these
62 * tracks and tracklets, and that histogram is stored in the output
65 * Furthermore, a number of diagnostics @f$\eta@f$ histograms are filled:
66 * - @f$\eta@f$ of all accepted tracks and tracklets
67 * - @f$\eta@f$ of accepted global tracks
68 * - @f$\eta@f$ of accepted ITS tracks
69 * - @f$\eta@f$ of accepted SPD tracklets
71 * At the end of the job, these histograms are normalized to the
72 * number of accepted events and bin width to provide an estimate of
73 * the @f$dN_{ch}/d\eta@f$
75 * Only minimum bias events with a @f$v_z@f$ within the defined cut
78 * @ingroup pwg2_forward_aod
80 class CentralMultTask : public AliAnalysisTaskSE
95 * @param name Name of task
96 * @param maxEta Set @f$\eta@f$ range
97 * @param maxVtx Set @f$v_z@f$ range
99 CentralMultTask(const char* name,
106 virtual ~CentralMultTask();
107 void SetUseTracklets(Bool_t use) { fUseTracklets = use; }
109 * Initialise on master - does nothing
112 virtual void Init() {}
114 * Create output objects.
116 * This is called once per slave process
118 virtual void UserCreateOutputObjects();
120 * Process a single event
122 * @param option Not used
124 virtual void UserExec(Option_t* option);
126 * Called at end of event processing.
128 * This is called once in the master
130 * @param option Not used
132 virtual void Terminate(Option_t* option);
135 * Check if this event is within cuts
138 * @param vz Vertex @f$ z@f$ coordinate
140 * @return true if this event is to be considered
142 UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
144 TH2D* fHist; // Per-event d2n/deta/dphi
145 TH1D* fAll; // Summed d2n/deta/dphi - all track(let)s
146 TH1D* fGlobal; // Summed d2n/deta/dphi - global tracks
147 TH1D* fITS; // Summed d2n/deta/dphi - ITS tracks
148 TH1D* fTracklets; // Summed d2n/deta/dphi - SPD tracklets
149 TH1D* fNEventsTr; // Number of triggered events per vertex bin
150 TH1D* fNEventsVtx;// Number of triggered+vertex events per vertex
151 TList* fOutput; // Output list
152 AliESDtrackCuts fQGlo; // Quality cut on ITS+TPC
153 AliESDtrackCuts fQITS; // Quality cut on ITS standalone (not ITS+TPC)
154 AliESDtrackCuts fDCAwSPD; // DCA for traks with SPD hits
155 AliESDtrackCuts fDCAwoSPD; // DCA for traks without SPD hits
156 AliESDtrackCuts fIsPrimary; // Primary particle
157 Bool_t fUseTracklets; // Whether to count tracklets or not
158 Bool_t fDebug; // Debug flag
160 ClassDef(CentralMultTask,1); // Determine multiplicity in central area
163 //====================================================================
169 #include <AliAnalysisManager.h>
170 #include <AliESDEvent.h>
171 #include <AliAODHandler.h>
172 #include <AliAODEvent.h>
173 #include <AliESDInputHandler.h>
174 #include <AliESDtrack.h>
175 #include <AliMultiplicity.h>
177 //____________________________________________________________________
178 inline CentralMultTask::CentralMultTask()
179 : AliAnalysisTaskSE(),
188 fUseTracklets(false),
192 //____________________________________________________________________
193 inline CentralMultTask::CentralMultTask(const char* /* name */,
196 : AliAnalysisTaskSE("Central"),
209 fHist = new TH2D("Central", "d^{2}N_{ch}/d#etad#phi in central region",
210 nBins, -maxEta, maxEta, 20, 0, 2*TMath::Pi());
211 fHist->SetDirectory(0);
212 fHist->SetXTitle("#eta");
213 fHist->SetYTitle("#phi [radians]");
214 fHist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
218 fAll = new TH1D("all", "Central region", nBins, -maxEta, maxEta);
219 fAll->SetDirectory(0);
220 fAll->SetXTitle("#eta");
221 fAll->SetYTitle("dN_{ch}/d#eta");
223 fAll->SetFillColor(kGray);
224 fAll->SetFillStyle(3001);
225 fAll->SetMarkerStyle(28);
226 fAll->SetMarkerColor(kGray);
229 fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
230 fGlobal->SetTitle("Global tracks");
231 fGlobal->SetDirectory(0);
233 fGlobal->SetFillColor(kRed+1);
234 fGlobal->SetFillStyle(3001);
235 fGlobal->SetMarkerStyle(28);
236 fGlobal->SetMarkerColor(kRed+1);
237 fGlobal->SetStats(0);
239 fITS = static_cast<TH1D*>(fAll->Clone("its"));
240 fITS->SetTitle("ITS tracks");
241 fITS->SetDirectory(0);
243 fITS->SetFillColor(kGreen+1);
244 fITS->SetFillStyle(3001);
245 fITS->SetMarkerStyle(28);
246 fITS->SetMarkerColor(kGreen+1);
249 fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
250 fTracklets->SetTitle("SPD tracklets");
251 fTracklets->SetDirectory(0);
253 fTracklets->SetFillColor(kBlue+1);
254 fTracklets->SetFillStyle(3001);
255 fTracklets->SetMarkerStyle(28);
256 fTracklets->SetMarkerColor(kBlue+1);
257 fTracklets->SetStats(0);
259 fNEventsTr = new TH1D("nEventsTr",
260 "Events per vertex bin", 10, -maxVtx, maxVtx);
261 fNEventsTr->SetXTitle("v_{z}");
262 fNEventsTr->SetYTitle("Events");
263 fNEventsTr->SetDirectory(0);
265 fNEventsVtx = new TH1D("nEventsVtx",
266 "Events per vertex bin", 10, -maxVtx, maxVtx);
267 fNEventsVtx->SetXTitle("v_{z}");
268 fNEventsVtx->SetYTitle("Events");
269 fNEventsVtx->SetDirectory(0);
271 // --- Global (ITS+TPC) track quality cuts
273 fQGlo.SetMinNClustersTPC(70);
274 fQGlo.SetMaxChi2PerClusterTPC(4);
275 fQGlo.SetAcceptKinkDaughters(kFALSE);
276 fQGlo.SetRequireTPCRefit(kTRUE);
278 fQGlo.SetRequireITSRefit(kTRUE);
279 fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
280 fQGlo.SetEtaRange(-maxEta, maxEta);
282 // --- ITS standalone track quality cuts
283 fQITS.SetRequireITSRefit(kTRUE);
284 fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
285 fQITS.SetEtaRange(-maxEta, maxEta);
287 // -- Distance-of-Closest-Approach cuts for tracks w/ITS hits
288 fDCAwSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
289 AliESDtrackCuts::kAny);
290 fDCAwSPD.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
291 fDCAwSPD.SetMaxDCAToVertexZ(0.5);
292 fDCAwSPD.SetEtaRange(-maxEta, maxEta);
294 // -- Distance-of-Closest-Approach cuts for tracks w/o ITS hits
295 fDCAwoSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
296 AliESDtrackCuts::kNone);
297 fDCAwoSPD.SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
298 fDCAwoSPD.SetMaxDCAToVertexZ(0.5);
299 fDCAwoSPD.SetEtaRange(-maxEta, maxEta);
301 // -- Primary track cut
302 // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
304 fIsPrimary.SetMinNClustersTPC(70);
305 fIsPrimary.SetMaxChi2PerClusterTPC(4);
306 fIsPrimary.SetAcceptKinkDaughters(kFALSE);
307 fIsPrimary.SetRequireTPCRefit(kTRUE);
308 fIsPrimary.SetRequireITSRefit(kTRUE);
309 fIsPrimary.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
310 AliESDtrackCuts::kAny);
312 fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
313 fIsPrimary.SetMaxDCAToVertexZ(2);
314 fIsPrimary.SetDCAToVertex2D(kFALSE);
315 fIsPrimary.SetRequireSigmaToVertex(kFALSE);
317 // Output slot #1 writes into a TH1 container
318 DefineOutput(1, TList::Class());
321 //____________________________________________________________________
322 inline CentralMultTask::~CentralMultTask()
324 if (fHist) delete fHist;
325 if (fAll) delete fAll;
326 if (fGlobal) delete fGlobal;
327 if (fITS) delete fITS;
328 if (fTracklets) delete fTracklets;
331 //________________________________________________________________________
333 CentralMultTask::UserCreateOutputObjects()
336 // Called once (on the worker node)
339 fOutput->SetName(GetName());
343 fOutput->Add(fGlobal);
345 fOutput->Add(fTracklets);
346 fOutput->Add(fNEventsTr);
347 fOutput->Add(fNEventsVtx);
349 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
351 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
352 if (!ah) AliFatal("No AOD output handler set in analysis manager");
354 ah->AddBranch("TH2D", &fHist);
356 // Post data for ALL output slots >0 here, to get at least an empty histogram
357 PostData(1, fOutput);
360 //____________________________________________________________________
362 CentralMultTask::CheckEvent(const AliESDEvent& esd, Double_t& vz)
364 // Do some fast cuts first
367 // Get the analysis manager - should always be there
368 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
370 AliWarning("No analysis manager defined!");
374 // Get the input handler - should always be there
375 AliInputEventHandler* ih =
376 static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
378 AliWarning("No input handler");
383 UInt_t mask = ih->IsEventSelected();
384 Bool_t isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
386 if (isMinBias) ret |= kValidTrigger;
388 // check for good reconstructed vertex
389 if (!(esd.GetPrimaryVertex()->GetStatus())) {
391 AliWarning("Primary vertex has bad status");
394 if (!(esd.GetPrimaryVertexSPD()->GetStatus())) {
396 AliWarning("Primary SPD vertex has bad status");
400 // if vertex is from spd vertexZ, require more stringent cut
401 if (esd.GetPrimaryVertex()->IsFromVertexerZ()) {
402 if (esd.GetPrimaryVertex()->GetDispersion() > 0.02 ||
403 esd.GetPrimaryVertex()->GetZRes() > 0.25) {
405 AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)",
406 esd.GetPrimaryVertex()->GetDispersion(),
407 esd.GetPrimaryVertex()->GetZRes()));
411 // One can add a cut on the vertex z position here
412 vz = esd.GetPrimaryVertex()->GetZ();
413 Double_t vl = fNEventsVtx->GetXaxis()->GetXmin();
414 Double_t vh = fNEventsVtx->GetXaxis()->GetXmax();
415 if (vz < vl || vz > vh) {
417 AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
425 //____________________________________________________________________
427 CentralMultTask::UserExec(Option_t *)
430 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
432 AliError("Cannot get the ESD event");
440 UShort_t eventFlags = CheckEvent(*esd, vz);
441 if (!(eventFlags & kValidTrigger))
442 return; // No trigger
444 fNEventsTr->Fill(vz);
445 if (eventFlags & kValidVertex)
446 fNEventsVtx->Fill(vz);
448 return; // No trigger or no vertex
452 // flags for secondary and rejected tracks
453 // set this bit in ESD tracks if it is rejected by a cut
454 const int kRejBit = BIT(15);
455 // set this bit in ESD tracks if it is secondary according to a cut
456 const int kSecBit = BIT(16);
459 Int_t nESDTracks = esd->GetNumberOfTracks();
460 // Loop over ESD tracks
461 for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
463 AliESDtrack* track = esd->GetTrack(i);
465 // if track is a secondary from a V0, flag as a secondary
466 if (track->IsOn(AliESDtrack::kMultInV0)) {
467 track->SetBit(kSecBit);
471 Double_t eta = track->Eta();
472 Double_t phi = track->Phi();
473 // check tracks with ITS part
474 if (track->IsOn(AliESDtrack::kITSin)) {
476 // Check ITS pure stand-alone - and if it is, reject it
477 if (track->IsOn(AliESDtrack::kITSpureSA)) {
478 track->SetBit(kRejBit);
482 // Check DCA - if not close enough reject it as secondary
483 if (!fDCAwSPD.AcceptTrack(track) &&
484 !fDCAwoSPD.AcceptTrack(track)) {
485 track->SetBit(kSecBit);
489 // Check if this is an ITS complementary - no TPC in
490 if (!track->IsOn(AliESDtrack::kTPCin)) {
491 if (fQITS.AcceptTrack(track)) { // Check ITS quality
496 track->SetBit(kRejBit);
498 else { // Not ITS SA or TPC+ITS
499 if (fQGlo.AcceptTrack(track)) { // Check ITS quality
504 track->SetBit(kRejBit);
506 if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
509 fHist->Fill(eta, phi);
513 // Get multiplicity from SPD tracklets
514 const AliMultiplicity* spdmult = esd->GetMultiplicity();
515 for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i){
516 // if counting tracks+tracklets,
517 // check if clusters were already used in tracks
519 spdmult->GetTrackletTrackIDs(i,0,id1,id2);
520 AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
521 AliESDtrack* tr2 = id2>=0 ? esd->GetTrack(id2) : 0;
523 if ((tr1 && tr1->TestBit(kSecBit)) || // Flagged as secondary
524 (tr2 && tr2->TestBit(kSecBit)) || // Flagged as secondary
525 (tr1 && !tr1->TestBit(kRejBit)) || // already accounted for
526 (tr2 && !tr2->TestBit(kRejBit))) // already accounted for
529 Double_t eta = spdmult->GetEta(i);
530 Double_t phi = spdmult->GetPhi(i);
532 // Increment d2n/detadphi from an SPD tracklet
534 fHist->Fill(eta, phi);
538 fTracklets->Fill(eta);
540 if (fDebug) AliInfo(Form("A total of %d tracks", total));
542 // NEW HISTO should be filled before this point, as PostData puts the
543 // information for this iteration of the UserExec in the container
544 PostData(1, fOutput);
548 //________________________________________________________________________
550 CentralMultTask::Terminate(Option_t *)
552 // Draw result to screen, or perform fitting, normalizations Called
553 // once at the end of the query
555 fOutput = dynamic_cast<TList*> (GetOutputData(1));
557 AliError("Could not retrieve TList fOutput");
561 TH1D* all = static_cast<TH1D*>(fOutput->FindObject("all"));
562 TH1D* global = static_cast<TH1D*>(fOutput->FindObject("global"));
563 TH1D* its = static_cast<TH1D*>(fOutput->FindObject("its"));
564 TH1D* tracklets = static_cast<TH1D*>(fOutput->FindObject("tracklets"));
565 TH1D* eventsTr = static_cast<TH1D*>(fOutput->FindObject("nEventsTr"));
566 TH1D* eventsVtx = static_cast<TH1D*>(fOutput->FindObject("nEventsVtx"));
568 Int_t nTriggers = eventsTr->GetEntries();
569 Int_t nVertex = eventsVtx->GetEntries();
570 if (nTriggers <= 0 || nVertex <= 0) {
571 AliWarning("No data in the events histogram!");
575 all ->Scale(1. / nTriggers, "width");
576 global ->Scale(1. / nTriggers, "width");
577 its ->Scale(1. / nTriggers, "width");
578 tracklets->Scale(1. / nTriggers, "width");
580 THStack* stack = new THStack("components", "Components");
581 if (fUseTracklets) stack->Add(tracklets);
588 //========================================================================
589 inline AliAnalysisTask*
590 AddTaskCentralTracks()
593 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
595 // Make our object. 2nd argumenent is absolute max Eta
596 // 3rd argument is absolute max Vz
597 CentralMultTask* task = new CentralMultTask("Global", 2, 10);
598 // if physics selection performed in UserExec(),
599 // this line should be commented
600 // task->SelectCollisionCandidates(AliVEvent::kMB);
603 // create containers for input/output
604 AliAnalysisDataContainer *output =
605 mgr->CreateContainer("Central", TList::Class(),
606 AliAnalysisManager::kOutputContainer,
607 AliAnalysisManager::GetCommonFileName());
609 // connect input/output
610 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
611 mgr->ConnectOutput(task, 1, output);
617 //________________________________________________________________________