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 pwglf_forward_scripts_tasks
12 #include <AliAnalysisTaskSE.h>
13 #include <AliESDtrackCuts.h>
19 * Task to determine the charged particle pseudo density from tracks
22 * This class determines
24 * \left.\frac{d^2N_{ch}}{d\eta d\phi}\right|_{central}
26 * from tracks and tracklets.
28 * First, global tracks are investigated. The requirements on global
30 * - @f$ n_{TPC clusters} \ge 70@f$
31 * - @f$ \chi^2_{TPC cluster} \le 4@f$
33 * - Re-fit of TPC tracks
34 * - Re-fit of ITS tracks
35 * - No requirement on SPD clusters
37 * Secondly, ITS stand-alone tracks are investigated. The
38 * requirements on ITS standalone tracks are
39 * - Re-fit of ITS tracks
40 * - No requirement on SPD clusters
42 * Tracks that does not meet these quality conditions are flagged as
45 * Both kinds of tracks (global and ITS standalone) have requirements
46 * on the distance of closest approach (DCA) to the interaction
47 * vertex @f$(v_x,v_y,v_z)@f$.
49 * For tracks with SPD clusters:
50 * - @f$ DCA_{xy} < 0.0182+0.0350/p_t^{1.01}@f$
51 * - @f$ DCA_{z} < 0.5@f$
53 * For tracks without SPD clusters
54 * - @f$ DCA_{xy} < 1.5(0.0182+0.0350/p_t^{1.01})@f$
55 * - @f$ DCA_{z} < 0.5@f$
57 * Tracks that does not meet these DCA requirements are flagged as
60 * Thirdly, the number of SPD tracklets are investigated. If the
61 * tracklet is associated with a track, and that track has already
62 * been used, then that tracklet is ignored
64 * An @f$(\eta,\phi)@f$ per-event histogram is then filled from these
65 * tracks and tracklets, and that histogram is stored in the output
68 * Furthermore, a number of diagnostics @f$\eta@f$ histograms are filled:
69 * - @f$\eta@f$ of all accepted tracks and tracklets
70 * - @f$\eta@f$ of accepted global tracks
71 * - @f$\eta@f$ of accepted ITS tracks
72 * - @f$\eta@f$ of accepted SPD tracklets
74 * At the end of the job, these histograms are normalized to the
75 * number of accepted events and bin width to provide an estimate of
76 * the @f$dN_{ch}/d\eta@f$
78 * Only minimum bias events with a @f$v_z@f$ within the defined cut
81 * @deprecated This class is deprecated
82 * @ingroup pwglf_forward_aod
84 class CentralMultTask : public AliAnalysisTaskSE
99 * @param name Name of task
100 * @param maxEta Set @f$\eta@f$ range
101 * @param maxVtx Set @f$v_z@f$ range
103 CentralMultTask(const char* name,
110 virtual ~CentralMultTask();
112 * Whether to use tracklets
114 * @param use IF true, use tracklets
116 void SetUseTracklets(Bool_t use) { fUseTracklets = use; }
118 * Initialise on master - does nothing
121 virtual void Init() {}
123 * Create output objects.
125 * This is called once per slave process
127 virtual void UserCreateOutputObjects();
129 * Process a single event
131 * @param option Not used
133 virtual void UserExec(Option_t* option);
135 * Called at end of event processing.
137 * This is called once in the master
139 * @param option Not used
141 virtual void Terminate(Option_t* option);
144 * Check if this event is within cuts
147 * @param vz Vertex @f$ z@f$ coordinate
149 * @return true if this event is to be considered
151 UShort_t CheckEvent(const AliESDEvent& esd, Double_t& vz);
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
169 ClassDef(CentralMultTask,1); // Determine multiplicity in central area
172 //====================================================================
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>
186 //____________________________________________________________________
187 inline CentralMultTask::CentralMultTask()
188 : AliAnalysisTaskSE(),
197 fUseTracklets(false),
201 //____________________________________________________________________
202 inline CentralMultTask::CentralMultTask(const char* /* name */,
205 : AliAnalysisTaskSE("Central"),
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");
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");
232 fAll->SetFillColor(kGray);
233 fAll->SetFillStyle(3001);
234 fAll->SetMarkerStyle(28);
235 fAll->SetMarkerColor(kGray);
238 fGlobal = static_cast<TH1D*>(fAll->Clone("global"));
239 fGlobal->SetTitle("Global tracks");
240 fGlobal->SetDirectory(0);
242 fGlobal->SetFillColor(kRed+1);
243 fGlobal->SetFillStyle(3001);
244 fGlobal->SetMarkerStyle(28);
245 fGlobal->SetMarkerColor(kRed+1);
246 fGlobal->SetStats(0);
248 fITS = static_cast<TH1D*>(fAll->Clone("its"));
249 fITS->SetTitle("ITS tracks");
250 fITS->SetDirectory(0);
252 fITS->SetFillColor(kGreen+1);
253 fITS->SetFillStyle(3001);
254 fITS->SetMarkerStyle(28);
255 fITS->SetMarkerColor(kGreen+1);
258 fTracklets = static_cast<TH1D*>(fAll->Clone("tracklets"));
259 fTracklets->SetTitle("SPD tracklets");
260 fTracklets->SetDirectory(0);
262 fTracklets->SetFillColor(kBlue+1);
263 fTracklets->SetFillStyle(3001);
264 fTracklets->SetMarkerStyle(28);
265 fTracklets->SetMarkerColor(kBlue+1);
266 fTracklets->SetStats(0);
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);
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);
280 // --- Global (ITS+TPC) track quality cuts
282 fQGlo.SetMinNClustersTPC(70);
283 fQGlo.SetMaxChi2PerClusterTPC(4);
284 fQGlo.SetAcceptKinkDaughters(kFALSE);
285 fQGlo.SetRequireTPCRefit(kTRUE);
287 fQGlo.SetRequireITSRefit(kTRUE);
288 fQGlo.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
289 fQGlo.SetEtaRange(-maxEta, maxEta);
291 // --- ITS standalone track quality cuts
292 fQITS.SetRequireITSRefit(kTRUE);
293 fQITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff);
294 fQITS.SetEtaRange(-maxEta, maxEta);
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);
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);
310 // -- Primary track cut
311 // https://twiki.cern.ch/twiki/bin/view/ALICE/SelectionOfPrimaryTracksForPpDataAnalysis
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);
321 fIsPrimary.SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
322 fIsPrimary.SetMaxDCAToVertexZ(2);
323 fIsPrimary.SetDCAToVertex2D(kFALSE);
324 fIsPrimary.SetRequireSigmaToVertex(kFALSE);
326 // Output slot #1 writes into a TH1 container
327 DefineOutput(1, TList::Class());
330 //____________________________________________________________________
331 inline CentralMultTask::~CentralMultTask()
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;
340 //________________________________________________________________________
342 CentralMultTask::UserCreateOutputObjects()
345 // Called once (on the worker node)
348 fOutput->SetName(GetName());
352 fOutput->Add(fGlobal);
354 fOutput->Add(fTracklets);
355 fOutput->Add(fNEventsTr);
356 fOutput->Add(fNEventsVtx);
358 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
360 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
361 if (!ah) AliFatal("No AOD output handler set in analysis manager");
363 ah->AddBranch("TH2D", &fHist);
365 // Post data for ALL output slots >0 here, to get at least an empty histogram
366 PostData(1, fOutput);
369 //____________________________________________________________________
371 CentralMultTask::CheckEvent(const AliESDEvent& esd, Double_t& vz)
373 // Do some fast cuts first
376 // Get the analysis manager - should always be there
377 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
379 AliWarning("No analysis manager defined!");
383 // Get the input handler - should always be there
384 AliInputEventHandler* ih =
385 static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
387 AliWarning("No input handler");
392 UInt_t mask = ih->IsEventSelected();
393 Bool_t isMinBias = (mask == AliVEvent::kMB) ? 1 : 0;
395 if (isMinBias) ret |= kValidTrigger;
397 // check for good reconstructed vertex
398 if (!(esd.GetPrimaryVertex()->GetStatus())) {
400 AliWarning("Primary vertex has bad status");
403 if (!(esd.GetPrimaryVertexSPD()->GetStatus())) {
405 AliWarning("Primary SPD vertex has bad status");
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) {
414 AliWarning(Form("Primary vertex dispersion=%f (0.02) zres=%f (0.05)",
415 esd.GetPrimaryVertex()->GetDispersion(),
416 esd.GetPrimaryVertex()->GetZRes()));
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) {
426 AliWarning(Form("Primary vertex vz=%f out of range [%f,%f]", vz, vl, vh));
434 //____________________________________________________________________
436 CentralMultTask::UserExec(Option_t *)
439 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
441 AliError("Cannot get the ESD event");
449 UShort_t eventFlags = CheckEvent(*esd, vz);
450 if (!(eventFlags & kValidTrigger))
451 return; // No trigger
453 fNEventsTr->Fill(vz);
454 if (eventFlags & kValidVertex)
455 fNEventsVtx->Fill(vz);
457 return; // No trigger or no vertex
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);
468 Int_t nESDTracks = esd->GetNumberOfTracks();
469 // Loop over ESD tracks
470 for(Int_t i = 0; i < nESDTracks; i++){ // flag the tracks
472 AliESDtrack* track = esd->GetTrack(i);
474 // if track is a secondary from a V0, flag as a secondary
475 if (track->IsOn(AliESDtrack::kMultInV0)) {
476 track->SetBit(kSecBit);
480 Double_t eta = track->Eta();
481 Double_t phi = track->Phi();
482 // check tracks with ITS part
483 if (track->IsOn(AliESDtrack::kITSin)) {
485 // Check ITS pure stand-alone - and if it is, reject it
486 if (track->IsOn(AliESDtrack::kITSpureSA)) {
487 track->SetBit(kRejBit);
491 // Check DCA - if not close enough reject it as secondary
492 if (!fDCAwSPD.AcceptTrack(track) &&
493 !fDCAwoSPD.AcceptTrack(track)) {
494 track->SetBit(kSecBit);
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
505 track->SetBit(kRejBit);
507 else { // Not ITS SA or TPC+ITS
508 if (fQGlo.AcceptTrack(track)) { // Check ITS quality
513 track->SetBit(kRejBit);
515 if (track->IsOn(kSecBit) || track->IsOn(kRejBit)) continue;
518 fHist->Fill(eta, phi);
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
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;
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
538 Double_t eta = spdmult->GetEta(i);
539 Double_t phi = spdmult->GetPhi(i);
541 // Increment d2n/detadphi from an SPD tracklet
543 fHist->Fill(eta, phi);
547 fTracklets->Fill(eta);
549 if (fDebug) AliInfo(Form("A total of %d tracks", total));
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);
557 //________________________________________________________________________
559 CentralMultTask::Terminate(Option_t *)
561 // Draw result to screen, or perform fitting, normalizations Called
562 // once at the end of the query
564 fOutput = dynamic_cast<TList*> (GetOutputData(1));
566 AliError("Could not retrieve TList fOutput");
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"));
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!");
584 all ->Scale(1. / nTriggers, "width");
585 global ->Scale(1. / nTriggers, "width");
586 its ->Scale(1. / nTriggers, "width");
587 tracklets->Scale(1. / nTriggers, "width");
589 THStack* stack = new THStack("components", "Components");
590 if (fUseTracklets) stack->Add(tracklets);
597 //========================================================================
599 * Add task to do central dN/deta from tracks/tracklets
601 * @deprecated This function is deprecated
602 * @return Newly created task
604 inline AliAnalysisTask*
605 AddTaskCentralTracks()
608 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
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);
618 // create containers for input/output
619 AliAnalysisDataContainer *output =
620 mgr->CreateContainer("Central", TList::Class(),
621 AliAnalysisManager::kOutputContainer,
622 AliAnalysisManager::GetCommonFileName());
624 // connect input/output
625 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
626 mgr->ConnectOutput(task, 1, output);
632 //________________________________________________________________________