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();
105 * Initialise on master - does nothing
108 virtual void Init() {}
110 * Create output objects.
112 * This is called once per slave process
114 virtual void UserCreateOutputObjects();
116 * Process a single event
118 * @param option Not used
120 virtual void UserExec(Option_t* option);
122 * Called at end of event processing..
124 * This is called once in the master
126 * @param option Not used
128 virtual void Terminate(Option_t* option);
131 * Check if this event is within cuts
135 * @return true if this event is to be considered
137 UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
139 TH2D* fHist; // Per-event d2n/deta/dphi
140 TH1D* fAll; // Summed d2n/deta/dphi - all track(let)s
141 TH1D* fGlobal; // Summed d2n/deta/dphi - global tracks
142 TH1D* fITS; // Summed d2n/deta/dphi - ITS tracks
143 TH1D* fTracklets; // Summed d2n/deta/dphi - SPD tracklets
144 TH1D* fNEventsTr; // Number of triggered events per vertex bin
145 TH1D* fNEventsVtx;// Number of triggered+vertex events per vertex
146 TList* fOutput; // Output list
147 AliESDtrackCuts fQGlo; // Quality cut on ITS+TPC
148 AliESDtrackCuts fQITS; // Quality cut on ITS standalone (not ITS+TPC)
149 AliESDtrackCuts fDCAwSPD; // DCA for traks with SPD hits
150 AliESDtrackCuts fDCAwoSPD; // DCA for traks without SPD hits
151 AliESDtrackCuts fIsPrimary; // Primary particle
152 Bool_t fDebug; // Debug flag
154 ClassDef(CentralMultTask,1); // Determine multiplicity in central area
157 //====================================================================
163 #include <AliAnalysisManager.h>
164 #include <AliESDEvent.h>
165 #include <AliAODHandler.h>
166 #include <AliAODEvent.h>
167 #include <AliESDInputHandler.h>
168 #include <AliESDtrack.h>
169 #include <AliMultiplicity.h>
171 //____________________________________________________________________
172 inline CentralMultTask::CentralMultTask()
173 : AliAnalysisTaskSE(),
185 //____________________________________________________________________
186 inline CentralMultTask::CentralMultTask(const char* name,
189 : AliAnalysisTaskSE(name),
200 fHist = new TH2D("Central", "d^{2}N_{ch}/d#etad#phi in central region",
201 80, -maxEta, maxEta, 20, 0, 2*TMath::Pi());
202 fHist->SetDirectory(0);
203 fHist->SetXTitle("#eta");
204 fHist->SetYTitle("#phi [radians]");
205 fHist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
209 fAll = new TH1D("all", "Central region", 80, -maxEta, maxEta);
210 fAll->SetDirectory(0);
211 fAll->SetXTitle("#eta");
212 fAll->SetYTitle("dN_{ch}/d#eta");
214 fAll->SetFillColor(kGray);
215 fAll->SetFillStyle(3001);
216 fAll->SetMarkerStyle(28);
217 fAll->SetMarkerColor(kGray);
220 fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
221 fGlobal->SetTitle("Global tracks");
222 fGlobal->SetDirectory(0);
224 fGlobal->SetFillColor(kRed+1);
225 fGlobal->SetFillStyle(3001);
226 fGlobal->SetMarkerStyle(28);
227 fGlobal->SetMarkerColor(kRed+1);
228 fGlobal->SetStats(0);
230 fITS = static_cast<TH1D*>(fAll->Clone("its"));
231 fITS->SetTitle("ITS tracks");
232 fITS->SetDirectory(0);
234 fITS->SetFillColor(kGreen+1);
235 fITS->SetFillStyle(3001);
236 fITS->SetMarkerStyle(28);
237 fITS->SetMarkerColor(kGreen+1);
240 fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
241 fTracklets->SetTitle("SPD tracklets");
242 fTracklets->SetDirectory(0);
244 fTracklets->SetFillColor(kBlue+1);
245 fTracklets->SetFillStyle(3001);
246 fTracklets->SetMarkerStyle(28);
247 fTracklets->SetMarkerColor(kBlue+1);
248 fTracklets->SetStats(0);
250 fNEventsTr = new TH1D("nEventsTr",
251 "Events per vertex bin", 10, -maxVtx, maxVtx);
252 fNEventsTr->SetXTitle("v_{z}");
253 fNEventsTr->SetYTitle("Events");
254 fNEventsTr->SetDirectory(0);
256 fNEventsVtx = new TH1D("nEventsVtx",
257 "Events per vertex bin", 10, -maxVtx, maxVtx);
258 fNEventsVtx->SetXTitle("v_{z}");
259 fNEventsVtx->SetYTitle("Events");
260 fNEventsVtx->SetDirectory(0);
262 // --- Global (ITS+TPC) track quality cuts
264 fQGlo.SetMinNClustersTPC(70);
265 fQGlo.SetMaxChi2PerClusterTPC(4);
266 fQGlo.SetAcceptKinkDaughters(kFALSE);
267 fQGlo.SetRequireTPCRefit(kTRUE);
269 fQGlo.SetRequireITSRefit(kTRUE);
270 fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
271 fQGlo.SetEtaRange(-maxEta, maxEta);
273 // --- ITS standalone track quality cuts
274 fQITS.SetRequireITSRefit(kTRUE);
275 fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
276 fQITS.SetEtaRange(-maxEta, maxEta);
278 // -- Distance-of-Closest-Approach cuts for tracks w/ITS hits
279 fDCAwSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
280 AliESDtrackCuts::kAny);
281 fDCAwSPD.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
282 fDCAwSPD.SetMaxDCAToVertexZ(0.5);
283 fDCAwSPD.SetEtaRange(-maxEta, maxEta);
285 // -- Distance-of-Closest-Approach cuts for tracks w/o ITS hits
286 fDCAwoSPD.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
287 AliESDtrackCuts::kNone);
288 fDCAwoSPD.SetMaxDCAToVertexXYPtDep("1.5*(0.0182+0.0350/pt^1.01)");
289 fDCAwoSPD.SetMaxDCAToVertexZ(0.5);
290 fDCAwoSPD.SetEtaRange(-maxEta, maxEta);
292 // -- Primary track cut
293 // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
295 fIsPrimary.SetMinNClustersTPC(70);
296 fIsPrimary.SetMaxChi2PerClusterTPC(4);
297 fIsPrimary.SetAcceptKinkDaughters(kFALSE);
298 fIsPrimary.SetRequireTPCRefit(kTRUE);
299 fIsPrimary.SetRequireITSRefit(kTRUE);
300 fIsPrimary.SetClusterRequirementITS(AliESDtrackCuts::kSPD,
301 AliESDtrackCuts::kAny);
303 fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
304 fIsPrimary.SetMaxDCAToVertexZ(2);
305 fIsPrimary.SetDCAToVertex2D(kFALSE);
306 fIsPrimary.SetRequireSigmaToVertex(kFALSE);
308 // Output slot #1 writes into a TH1 container
309 DefineOutput(1, TList::Class());
312 //____________________________________________________________________
313 inline CentralMultTask::~CentralMultTask()
315 if (fHist) delete fHist;
316 if (fAll) delete fAll;
317 if (fGlobal) delete fGlobal;
318 if (fITS) delete fITS;
319 if (fTracklets) delete fTracklets;
322 //________________________________________________________________________
324 CentralMultTask::UserCreateOutputObjects()
327 // Called once (on the worker node)
330 fOutput->SetName(GetName());
334 fOutput->Add(fGlobal);
336 fOutput->Add(fTracklets);
337 fOutput->Add(fNEventsTr);
338 fOutput->Add(fNEventsVtx);
340 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
342 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
343 if (!ah) AliFatal("No AOD output handler set in analysis manager");
345 ah->AddBranch("TH2D", &fHist);
347 // Post data for ALL output slots >0 here, to get at least an empty histogram
348 PostData(1, fOutput);
351 //____________________________________________________________________
353 CentralMultTask::CheckEvent(const AliESDEvent& esd, Double_t& vz)
355 // Do some fast cuts first
358 // Get the analysis manager - should always be there
359 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
361 AliWarning("No analysis manager defined!");
365 // Get the input handler - should always be there
366 AliInputEventHandler* ih =
367 static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
369 AliWarning("No input handler");
374 UInt_t mask = ih->IsEventSelected();
375 Bool_t isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
377 if (isMinBias) ret |= kValidTrigger;
379 // check for good reconstructed vertex
380 if (!(esd.GetPrimaryVertex()->GetStatus())) {
382 AliWarning("Primary vertex has bad status");
385 if (!(esd.GetPrimaryVertexSPD()->GetStatus())) {
387 AliWarning("Primary SPD vertex has bad status");
391 // if vertex is from spd vertexZ, require more stringent cut
392 if (esd.GetPrimaryVertex()->IsFromVertexerZ()) {
393 if (esd.GetPrimaryVertex()->GetDispersion() > 0.02 ||
394 esd.GetPrimaryVertex()->GetZRes() > 0.25) {
396 AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)",
397 esd.GetPrimaryVertex()->GetDispersion(),
398 esd.GetPrimaryVertex()->GetZRes()));
402 // One can add a cut on the vertex z position here
403 vz = esd.GetPrimaryVertex()->GetZ();
404 Double_t vl = fNEventsVtx->GetXaxis()->GetXmin();
405 Double_t vh = fNEventsVtx->GetXaxis()->GetXmax();
406 if (vz < vl || vz > vh) {
408 AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
416 //____________________________________________________________________
418 CentralMultTask::UserExec(Option_t *)
421 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
423 AliError("Cannot get the ESD event");
431 UShort_t eventFlags = CheckEvent(*esd, vz);
432 if (!(eventFlags & kValidTrigger))
433 return; // No trigger
435 fNEventsTr->Fill(vz);
436 if (eventFlags & kValidVertex)
437 fNEventsVtx->Fill(vz);
439 return; // No trigger or no vertex
443 // flags for secondary and rejected tracks
444 // set this bit in ESD tracks if it is rejected by a cut
445 const int kRejBit = BIT(15);
446 // set this bit in ESD tracks if it is secondary according to a cut
447 const int kSecBit = BIT(16);
450 Int_t nESDTracks = esd->GetNumberOfTracks();
451 // Loop over ESD tracks
452 for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
454 AliESDtrack* track = esd->GetTrack(i);
456 // if track is a secondary from a V0, flag as a secondary
457 if (track->IsOn(AliESDtrack::kMultInV0)) {
458 track->SetBit(kSecBit);
462 Double_t eta = track->Eta();
463 Double_t phi = track->Phi();
464 // check tracks with ITS part
465 if (track->IsOn(AliESDtrack::kITSin) &&
466 !track->IsOn(AliESDtrack::kITSpureSA)) {
467 // track has ITS part but is not an ITS_SA -> TPC+ITS
468 if (track->IsOn(AliESDtrack::kTPCin)) {
469 // Global track, has ITS and TPC contributions
470 if (fQGlo.AcceptTrack(track)) { // good TPCITS track
471 if (fDCAwSPD.AcceptTrack(track) ||
472 fDCAwoSPD.AcceptTrack(track)) {
478 // large DCA -> secondary, don't count either track not
479 // associated tracklet
480 track->SetBit(kSecBit);
483 // bad quality, don't count the track, but may count
484 // tracklet if associated
485 track->SetBit(kRejBit);
486 } // if (track->IsOn(AliESDtrack::kTPCin))
488 else if (fQITS.AcceptTrack(track)) {
489 // good ITS complimentary track
490 if (fDCAwSPD.AcceptTrack(track) ||
491 fDCAwoSPD.AcceptTrack(track)) {
497 // large DCA -> secondary, don't count either track not
498 // associated tracklet
499 track->SetBit(kSecBit);
500 } // else if (fQITS.AcceptTrack(track))
502 // bad quality, don't count the track, but may count tracklet
504 track->SetBit(kRejBit);
505 if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
508 fHist->Fill(eta, phi);
509 } /* track->IsOn(AliESDtrack::kITSin) &&
510 * !track->IsOn(AliESDtrack::kITSpureSA) */
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)) ||
524 (tr2 && tr2->TestBit(kSecBit)))
525 // reject as secondary
527 if ((tr1 && !tr1->TestBit(kRejBit)) ||
528 (tr2 && !tr2->TestBit(kRejBit)))
529 // already accounted for
532 Double_t eta = spdmult->GetEta(i);
533 Double_t phi = spdmult->GetPhi(i);
535 // Incremetn d2n/detadphi from an SPD tracklet
536 fHist->Fill(eta, phi);
537 fTracklets->Fill(eta);
541 AliInfo(Form("A total of %d tracks", total));
543 // NEW HISTO should be filled before this point, as PostData puts the
544 // information for this iteration of the UserExec in the container
545 PostData(1, fOutput);
549 //________________________________________________________________________
551 CentralMultTask::Terminate(Option_t *)
553 // Draw result to screen, or perform fitting, normalizations Called
554 // once at the end of the query
556 fOutput = dynamic_cast<TList*> (GetOutputData(1));
558 AliError("Could not retrieve TList fOutput");
562 TH1D* all = static_cast<TH1D*>(fOutput->FindObject("all"));
563 TH1D* global = static_cast<TH1D*>(fOutput->FindObject("global"));
564 TH1D* its = static_cast<TH1D*>(fOutput->FindObject("its"));
565 TH1D* tracklets = static_cast<TH1D*>(fOutput->FindObject("tracklets"));
566 TH1D* eventsTr = static_cast<TH1D*>(fOutput->FindObject("nEventsTr"));
567 TH1D* eventsVtx = static_cast<TH1D*>(fOutput->FindObject("nEventsVtx"));
569 Int_t nTriggers = eventsTr->GetEntries();
570 Int_t nVertex = eventsVtx->GetEntries();
571 if (nTriggers <= 0 || nVertex <= 0) {
572 AliWarning("No data in the events histogram!");
576 all ->Scale(1. / nTriggers, "width");
577 global ->Scale(1. / nTriggers, "width");
578 its ->Scale(1. / nTriggers, "width");
579 tracklets->Scale(1. / nTriggers, "width");
581 THStack* stack = new THStack("components", "Components");
582 stack->Add(tracklets);
589 //========================================================================
590 inline AliAnalysisTask*
594 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
596 // Make our object. 2nd argumenent is absolute max Eta
597 // 3rd argument is absolute max Vz
598 CentralMultTask* task = new CentralMultTask("Global", 2, 10);
599 // if physics selection performed in UserExec(),
600 // this line should be commented
601 // task->SelectCollisionCandidates(AliVEvent::kMB);
604 // create containers for input/output
605 AliAnalysisDataContainer *output =
606 mgr->CreateContainer("Central", TList::Class(),
607 AliAnalysisManager::kOutputContainer,
608 AliAnalysisManager::GetCommonFileName());
610 // connect input/output
611 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
612 mgr->ConnectOutput(task, 1, output);
618 //________________________________________________________________________