2 * @file AddCentralMultTask.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
11 #include <AliAnalysisTaskSE.h>
12 #include <AliESDtrackCuts.h>
18 * Task to determine the
20 * \left.\frac{d^2N_{ch}}{d\eta d\phi}\right|_{central}
22 * from tracks and tracklets.
24 * First, global tracks are investigated. The requirements on global
26 * - @f$ n_{TPC clusters} \ge 70@f$
27 * - @f$ \chi^2_{TPC cluster} \le 4@f$
29 * - Re-fit of TPC tracks
30 * - Re-fit of ITS tracks
31 * - No requirement on SPD clusters
33 * Secondly, ITS stand-alone tracks are investigated. The
34 * requirements on ITS standalone tracks are
35 * - Re-fit of ITS tracks
36 * - No requirement on SPD clusters
38 * Tracks that does not meet these quality conditions are flagged as
41 * Both kinds of tracks (global and ITS standalone) have requirements
42 * on the distance of closest approach (DCA) to the interaction
43 * vertex @f$(v_x,v_y,v_z)@f$.
45 * For tracks with SPD clusters:
46 * - @f$ DCA_{xy} < 0.0182+0.0350/p_t^{1.01}@f$
47 * - @f$ DCA_{z} < 0.5@f$
49 * For tracks without SPD clusters
50 * - @f$ DCA_{xy} < 1.5(0.0182+0.0350/p_t^{1.01})@f$
51 * - @f$ DCA_{z} < 0.5@f$
53 * Tracks that does not meet these DCA requirements are flagged as
56 * Thirdly, the number of SPD tracklets are investigated. If the
57 * tracklet is associated with a track, and that track has already
58 * been used, then that tracklet is ignored
60 * An @f$(\eta,\phi)@f$ per-event histogram is then filled from these
61 * tracks and tracklets, and that histogram is stored in the output
64 * Furthermore, a number of diagnostics @f$\eta@f$ histograms are filled:
65 * - @f$\eta@f$ of all accepted tracks and tracklets
66 * - @f$\eta@f$ of accepted global tracks
67 * - @f$\eta@f$ of accepted ITS tracks
68 * - @f$\eta@f$ of accepted SPD tracklets
70 * At the end of the job, these histograms are normalized to the
71 * number of accepted events and bin width to provide an estimate of
72 * the @f$dN_{ch}/d\eta@f$
74 * Only minimum bias events with a @f$v_z@f$ within the defined cut
77 class CentralMultTask : public AliAnalysisTaskSE
92 * @param name Name of task
93 * @param maxEta Set @f$\eta@f$ range
94 * @param maxVtx Set @f$v_z@f$ range
96 CentralMultTask(const char* name,
103 virtual ~CentralMultTask();
104 void SetUseTracklets(Bool_t use) { fUseTracklets = use; }
106 * Initialise on master - does nothing
109 virtual void Init() {}
111 * Create output objects.
113 * This is called once per slave process
115 virtual void UserCreateOutputObjects();
117 * Process a single event
119 * @param option Not used
121 virtual void UserExec(Option_t* option);
123 * Called at end of event processing..
125 * This is called once in the master
127 * @param option Not used
129 virtual void Terminate(Option_t* option);
132 * Check if this event is within cuts
136 * @return true if this event is to be considered
138 UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
140 TH2D* fHist; // Per-event d2n/deta/dphi
141 TH1D* fAll; // Summed d2n/deta/dphi - all track(let)s
142 TH1D* fGlobal; // Summed d2n/deta/dphi - global tracks
143 TH1D* fITS; // Summed d2n/deta/dphi - ITS tracks
144 TH1D* fTracklets; // Summed d2n/deta/dphi - SPD tracklets
145 TH1D* fNEventsTr; // Number of triggered events per vertex bin
146 TH1D* fNEventsVtx;// Number of triggered+vertex events per vertex
147 TList* fOutput; // Output list
148 AliESDtrackCuts fQGlo; // Quality cut on ITS+TPC
149 AliESDtrackCuts fQITS; // Quality cut on ITS standalone (not ITS+TPC)
150 AliESDtrackCuts fDCAwSPD; // DCA for traks with SPD hits
151 AliESDtrackCuts fDCAwoSPD; // DCA for traks without SPD hits
152 AliESDtrackCuts fIsPrimary; // Primary particle
153 Bool_t fUseTracklets; // Whether to count tracklets or not
154 Bool_t fDebug; // Debug flag
156 ClassDef(CentralMultTask,1); // Determine multiplicity in central area
159 //====================================================================
165 #include <AliAnalysisManager.h>
166 #include <AliESDEvent.h>
167 #include <AliAODHandler.h>
168 #include <AliAODEvent.h>
169 #include <AliESDInputHandler.h>
170 #include <AliESDtrack.h>
171 #include <AliMultiplicity.h>
173 //____________________________________________________________________
174 inline CentralMultTask::CentralMultTask()
175 : AliAnalysisTaskSE(),
184 fUseTracklets(false),
188 //____________________________________________________________________
189 inline CentralMultTask::CentralMultTask(const char* /* name */,
192 : AliAnalysisTaskSE("Central"),
205 fHist = new TH2D("Central", "d^{2}N_{ch}/d#etad#phi in central region",
206 nBins, -maxEta, maxEta, 20, 0, 2*TMath::Pi());
207 fHist->SetDirectory(0);
208 fHist->SetXTitle("#eta");
209 fHist->SetYTitle("#phi [radians]");
210 fHist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
214 fAll = new TH1D("all", "Central region", nBins, -maxEta, maxEta);
215 fAll->SetDirectory(0);
216 fAll->SetXTitle("#eta");
217 fAll->SetYTitle("dN_{ch}/d#eta");
219 fAll->SetFillColor(kGray);
220 fAll->SetFillStyle(3001);
221 fAll->SetMarkerStyle(28);
222 fAll->SetMarkerColor(kGray);
225 fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
226 fGlobal->SetTitle("Global tracks");
227 fGlobal->SetDirectory(0);
229 fGlobal->SetFillColor(kRed+1);
230 fGlobal->SetFillStyle(3001);
231 fGlobal->SetMarkerStyle(28);
232 fGlobal->SetMarkerColor(kRed+1);
233 fGlobal->SetStats(0);
235 fITS = static_cast<TH1D*>(fAll->Clone("its"));
236 fITS->SetTitle("ITS tracks");
237 fITS->SetDirectory(0);
239 fITS->SetFillColor(kGreen+1);
240 fITS->SetFillStyle(3001);
241 fITS->SetMarkerStyle(28);
242 fITS->SetMarkerColor(kGreen+1);
245 fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
246 fTracklets->SetTitle("SPD tracklets");
247 fTracklets->SetDirectory(0);
249 fTracklets->SetFillColor(kBlue+1);
250 fTracklets->SetFillStyle(3001);
251 fTracklets->SetMarkerStyle(28);
252 fTracklets->SetMarkerColor(kBlue+1);
253 fTracklets->SetStats(0);
255 fNEventsTr = new TH1D("nEventsTr",
256 "Events per vertex bin", 10, -maxVtx, maxVtx);
257 fNEventsTr->SetXTitle("v_{z}");
258 fNEventsTr->SetYTitle("Events");
259 fNEventsTr->SetDirectory(0);
261 fNEventsVtx = new TH1D("nEventsVtx",
262 "Events per vertex bin", 10, -maxVtx, maxVtx);
263 fNEventsVtx->SetXTitle("v_{z}");
264 fNEventsVtx->SetYTitle("Events");
265 fNEventsVtx->SetDirectory(0);
267 // --- Global (ITS+TPC) track quality cuts
269 fQGlo.SetMinNClustersTPC(70);
270 fQGlo.SetMaxChi2PerClusterTPC(4);
271 fQGlo.SetAcceptKinkDaughters(kFALSE);
272 fQGlo.SetRequireTPCRefit(kTRUE);
274 fQGlo.SetRequireITSRefit(kTRUE);
275 fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
276 fQGlo.SetEtaRange(-maxEta, maxEta);
278 // --- ITS standalone track quality cuts
279 fQITS.SetRequireITSRefit(kTRUE);
280 fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
281 fQITS.SetEtaRange(-maxEta, maxEta);
283 // -- Distance-of-Closest-Approach cuts for tracks w/ITS hits
284 fDCAwSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
285 AliESDtrackCuts::kAny);
286 fDCAwSPD.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
287 fDCAwSPD.SetMaxDCAToVertexZ(0.5);
288 fDCAwSPD.SetEtaRange(-maxEta, maxEta);
290 // -- Distance-of-Closest-Approach cuts for tracks w/o ITS hits
291 fDCAwoSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
292 AliESDtrackCuts::kNone);
293 fDCAwoSPD.SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
294 fDCAwoSPD.SetMaxDCAToVertexZ(0.5);
295 fDCAwoSPD.SetEtaRange(-maxEta, maxEta);
297 // -- Primary track cut
298 // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
300 fIsPrimary.SetMinNClustersTPC(70);
301 fIsPrimary.SetMaxChi2PerClusterTPC(4);
302 fIsPrimary.SetAcceptKinkDaughters(kFALSE);
303 fIsPrimary.SetRequireTPCRefit(kTRUE);
304 fIsPrimary.SetRequireITSRefit(kTRUE);
305 fIsPrimary.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
306 AliESDtrackCuts::kAny);
308 fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
309 fIsPrimary.SetMaxDCAToVertexZ(2);
310 fIsPrimary.SetDCAToVertex2D(kFALSE);
311 fIsPrimary.SetRequireSigmaToVertex(kFALSE);
313 // Output slot #1 writes into a TH1 container
314 DefineOutput(1, TList::Class());
317 //____________________________________________________________________
318 inline CentralMultTask::~CentralMultTask()
320 if (fHist) delete fHist;
321 if (fAll) delete fAll;
322 if (fGlobal) delete fGlobal;
323 if (fITS) delete fITS;
324 if (fTracklets) delete fTracklets;
327 //________________________________________________________________________
329 CentralMultTask::UserCreateOutputObjects()
332 // Called once (on the worker node)
335 fOutput->SetName(GetName());
339 fOutput->Add(fGlobal);
341 fOutput->Add(fTracklets);
342 fOutput->Add(fNEventsTr);
343 fOutput->Add(fNEventsVtx);
345 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
347 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
348 if (!ah) AliFatal("No AOD output handler set in analysis manager");
350 ah->AddBranch("TH2D", &fHist);
352 // Post data for ALL output slots >0 here, to get at least an empty histogram
353 PostData(1, fOutput);
356 //____________________________________________________________________
358 CentralMultTask::CheckEvent(const AliESDEvent& esd, Double_t& vz)
360 // Do some fast cuts first
363 // Get the analysis manager - should always be there
364 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
366 AliWarning("No analysis manager defined!");
370 // Get the input handler - should always be there
371 AliInputEventHandler* ih =
372 static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
374 AliWarning("No input handler");
379 UInt_t mask = ih->IsEventSelected();
380 Bool_t isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
382 if (isMinBias) ret |= kValidTrigger;
384 // check for good reconstructed vertex
385 if (!(esd.GetPrimaryVertex()->GetStatus())) {
387 AliWarning("Primary vertex has bad status");
390 if (!(esd.GetPrimaryVertexSPD()->GetStatus())) {
392 AliWarning("Primary SPD vertex has bad status");
396 // if vertex is from spd vertexZ, require more stringent cut
397 if (esd.GetPrimaryVertex()->IsFromVertexerZ()) {
398 if (esd.GetPrimaryVertex()->GetDispersion() > 0.02 ||
399 esd.GetPrimaryVertex()->GetZRes() > 0.25) {
401 AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)",
402 esd.GetPrimaryVertex()->GetDispersion(),
403 esd.GetPrimaryVertex()->GetZRes()));
407 // One can add a cut on the vertex z position here
408 vz = esd.GetPrimaryVertex()->GetZ();
409 Double_t vl = fNEventsVtx->GetXaxis()->GetXmin();
410 Double_t vh = fNEventsVtx->GetXaxis()->GetXmax();
411 if (vz < vl || vz > vh) {
413 AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
421 //____________________________________________________________________
423 CentralMultTask::UserExec(Option_t *)
426 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
428 AliError("Cannot get the ESD event");
436 UShort_t eventFlags = CheckEvent(*esd, vz);
437 if (!(eventFlags & kValidTrigger))
438 return; // No trigger
440 fNEventsTr->Fill(vz);
441 if (eventFlags & kValidVertex)
442 fNEventsVtx->Fill(vz);
444 return; // No trigger or no vertex
448 // flags for secondary and rejected tracks
449 // set this bit in ESD tracks if it is rejected by a cut
450 const int kRejBit = BIT(15);
451 // set this bit in ESD tracks if it is secondary according to a cut
452 const int kSecBit = BIT(16);
455 Int_t nESDTracks = esd->GetNumberOfTracks();
456 // Loop over ESD tracks
457 for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
459 AliESDtrack* track = esd->GetTrack(i);
461 // if track is a secondary from a V0, flag as a secondary
462 if (track->IsOn(AliESDtrack::kMultInV0)) {
463 track->SetBit(kSecBit);
467 Double_t eta = track->Eta();
468 Double_t phi = track->Phi();
469 // check tracks with ITS part
470 if (track->IsOn(AliESDtrack::kITSin)) {
472 // Check ITS pure stand-alone - and if it is, reject it
473 if (track->IsOn(AliESDtrack::kITSpureSA)) {
474 track->SetBit(kRejBit);
478 // Check DCA - if not close enough reject it as secondary
479 if (!fDCAwSPD.AcceptTrack(track) &&
480 !fDCAwoSPD.AcceptTrack(track)) {
481 track->SetBit(kSecBit);
485 // Check if this is an ITS complementary - no TPC in
486 if (!track->IsOn(AliESDtrack::kTPCin)) {
487 if (fQITS.AcceptTrack(track)) { // Check ITS quality
492 track->SetBit(kRejBit);
494 else { // Not ITS SA or TPC+ITS
495 if (fQGlo.AcceptTrack(track)) { // Check ITS quality
500 track->SetBit(kRejBit);
502 if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
505 fHist->Fill(eta, phi);
509 // Get multiplicity from SPD tracklets
510 const AliMultiplicity* spdmult = esd->GetMultiplicity();
511 for (Int_t i=0; i<spdmult->GetNumberOfTracklets(); ++i){
512 // if counting tracks+tracklets,
513 // check if clusters were already used in tracks
515 spdmult->GetTrackletTrackIDs(i,0,id1,id2);
516 AliESDtrack* tr1 = id1>=0 ? esd->GetTrack(id1) : 0;
517 AliESDtrack* tr2 = id2>=0 ? esd->GetTrack(id2) : 0;
519 if ((tr1 && tr1->TestBit(kSecBit)) || // Flagged as secondary
520 (tr2 && tr2->TestBit(kSecBit)) || // Flagged as secondary
521 (tr1 && !tr1->TestBit(kRejBit)) || // already accounted for
522 (tr2 && !tr2->TestBit(kRejBit))) // already accounted for
525 Double_t eta = spdmult->GetEta(i);
526 Double_t phi = spdmult->GetPhi(i);
528 // Increment d2n/detadphi from an SPD tracklet
530 fHist->Fill(eta, phi);
534 fTracklets->Fill(eta);
536 if (fDebug) AliInfo(Form("A total of %d tracks", total));
538 // NEW HISTO should be filled before this point, as PostData puts the
539 // information for this iteration of the UserExec in the container
540 PostData(1, fOutput);
544 //________________________________________________________________________
546 CentralMultTask::Terminate(Option_t *)
548 // Draw result to screen, or perform fitting, normalizations Called
549 // once at the end of the query
551 fOutput = dynamic_cast<TList*> (GetOutputData(1));
553 AliError("Could not retrieve TList fOutput");
557 TH1D* all = static_cast<TH1D*>(fOutput->FindObject("all"));
558 TH1D* global = static_cast<TH1D*>(fOutput->FindObject("global"));
559 TH1D* its = static_cast<TH1D*>(fOutput->FindObject("its"));
560 TH1D* tracklets = static_cast<TH1D*>(fOutput->FindObject("tracklets"));
561 TH1D* eventsTr = static_cast<TH1D*>(fOutput->FindObject("nEventsTr"));
562 TH1D* eventsVtx = static_cast<TH1D*>(fOutput->FindObject("nEventsVtx"));
564 Int_t nTriggers = eventsTr->GetEntries();
565 Int_t nVertex = eventsVtx->GetEntries();
566 if (nTriggers <= 0 || nVertex <= 0) {
567 AliWarning("No data in the events histogram!");
571 all ->Scale(1. / nTriggers, "width");
572 global ->Scale(1. / nTriggers, "width");
573 its ->Scale(1. / nTriggers, "width");
574 tracklets->Scale(1. / nTriggers, "width");
576 THStack* stack = new THStack("components", "Components");
577 if (fUseTracklets) stack->Add(tracklets);
584 //========================================================================
585 inline AliAnalysisTask*
589 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
591 // Make our object. 2nd argumenent is absolute max Eta
592 // 3rd argument is absolute max Vz
593 CentralMultTask* task = new CentralMultTask("Global", 2, 10);
594 // if physics selection performed in UserExec(),
595 // this line should be commented
596 // task->SelectCollisionCandidates(AliVEvent::kMB);
599 // create containers for input/output
600 AliAnalysisDataContainer *output =
601 mgr->CreateContainer("Central", TList::Class(),
602 AliAnalysisManager::kOutputContainer,
603 AliAnalysisManager::GetCommonFileName());
605 // connect input/output
606 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
607 mgr->ConnectOutput(task, 1, output);
613 //________________________________________________________________________