1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // Analysis task for Kinematic filtering
18 // Fill AOD tracks from Kinematic stack
20 // Code taken from macro CreateAODfromKineTree.C by Markus Oldenburg
26 #include "AliAnalysisTaskKineFilter.h"
27 #include "AliAnalysisManager.h"
28 #include "AliAnalysisFilter.h"
29 #include "AliHeader.h"
31 #include "TParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliAODEvent.h"
34 #include "AliAODHeader.h"
35 #include "AliAODVertex.h"
36 #include "AliAODTrack.h"
40 ClassImp(AliAnalysisTaskKineFilter)
42 ////////////////////////////////////////////////////////////////////////
44 //____________________________________________________________________
45 AliAnalysisTaskKineFilter::AliAnalysisTaskKineFilter():
48 // Default constructor
51 //____________________________________________________________________
52 AliAnalysisTaskKineFilter::AliAnalysisTaskKineFilter(const char* name):
53 AliAnalysisTaskSE(name),
56 // Default constructor
57 DefineInput (0, TChain::Class());
58 DefineOutput(0, TTree::Class());
61 //____________________________________________________________________
62 AliAnalysisTaskKineFilter::AliAnalysisTaskKineFilter(const AliAnalysisTaskKineFilter& obj):
63 AliAnalysisTaskSE(obj),
67 fTrackFilter = obj.fTrackFilter;
69 //____________________________________________________________________
70 AliAnalysisTaskKineFilter::~AliAnalysisTaskKineFilter()
72 // if( fTrackFilter ) delete fTrackFilter;
76 //____________________________________________________________________
77 AliAnalysisTaskKineFilter& AliAnalysisTaskKineFilter::operator=(const AliAnalysisTaskKineFilter& other)
80 AliAnalysisTaskSE::operator=(other);
81 fTrackFilter = other.fTrackFilter;
85 //____________________________________________________________________
86 void AliAnalysisTaskKineFilter::UserCreateOutputObjects()
88 // Create the output container
90 OutputTree()->GetUserInfo()->Add(fTrackFilter);
94 //____________________________________________________________________
95 void AliAnalysisTaskKineFilter::UserExec(Option_t */*option*/)
97 // Execute analysis for current event
100 // Fill AOD tracks from Kinematic stack
103 AliAODEvent* aod = AODEvent();
105 AliWarning("No Output Handler connected, doing nothing !") ;
109 // aod->CreateStdContent();
111 AliStack* stack = MCEvent()->Stack();
112 Int_t nTracks = stack->GetNtrack();
113 Int_t nPrims = stack->GetNprimary();
116 AliAODVertex *primary = NULL;
124 // Access to the header
125 AliAODHeader *header = aod->GetHeader();
127 Double_t emEnergy[2] = {-999., -999.};
130 *header = AliAODHeader(MCEvent()->Header()->GetRun(),
138 -999., // muon mag. field
140 -999., // ZDCN1Energy
141 -999., // ZDCP1Energy
142 -999., // ZDCN2Energy
143 -999., // ZDCP2Energy
144 emEnergy, // emEnergy
150 // Access to the AOD container of vertices
151 TClonesArray &vertices = *(aod->GetVertices());
153 // Access to the AOD container of tracks
154 TClonesArray &tracks = *(aod->GetTracks());
156 aod->ResetStd(nTracks, 1);
159 for (Int_t iTrack = 0; iTrack < nPrims; ++iTrack) {
161 TParticle *part = stack->Particle(iTrack);
164 // add primary vertex
165 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
166 primary = new(vertices[jVertices++])
167 AliAODVertex(x, NULL, -999., NULL, -1, AliAODVertex::kPrimary);
173 UInt_t selectInfo = 0;
175 selectInfo = fTrackFilter->IsSelected(part);
176 if (!selectInfo) continue;
179 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
180 p[0] = part->Px(); p[1] = part->Py(); p[2] = part->Pz();
182 // add primary tracks
183 primary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID,
191 0, // no ITSClusterMap
194 kFALSE, // no fit performed
195 kFALSE, // no fit preformed
196 AliAODTrack::kPrimary,
199 AliAODTrack* currTrack = (AliAODTrack*)tracks.Last();
200 SetChargeAndPID(part->GetPdgCode(), currTrack);
201 if (currTrack->Charge() != -99) {
202 if (currTrack->Charge() > 0) {
204 } else if (currTrack->Charge() < 0) {
209 LoopOverSecondaries(part, jTracks, jVertices, nPos, nNeg);
211 } // end of track loop
213 header->SetRefMultiplicityPos(nPos);
214 header->SetRefMultiplicityNeg(nNeg);
218 AliInfo(Form("primaries: %d secondaries: %d (pos: %d neg: %d), vertices: %d",
219 nPrimsAdd, tracks.GetEntriesFast()-nPrimsAdd, nPos, nNeg, vertices.GetEntriesFast() ) );
224 //____________________________________________________________________
225 Int_t AliAnalysisTaskKineFilter::LoopOverSecondaries(TParticle *mother, Int_t &jTracks, Int_t &jVertices, Int_t &nPos, Int_t &nNeg )
228 if (mother->GetNDaughters() > 0) {
230 AliStack* stack = MCEvent()->Stack();
232 TClonesArray &vertices = *(AODEvent()->GetVertices());
233 TClonesArray &tracks = *(AODEvent()->GetTracks());
236 AliAODVertex* secondary = NULL;
238 for (Int_t iDaughter = mother->GetFirstDaughter(); iDaughter <= mother->GetLastDaughter(); iDaughter++) {
239 TParticle *part = stack->Particle(iDaughter);
240 // only final particles
249 if (iDaughter == mother->GetFirstDaughter()) {
250 // add secondary vertex
251 secondary = new(vertices[jVertices++])
252 AliAODVertex(x, NULL, -999., tracks.Last(), iDaughter, AliAODVertex::kUndef);
253 SetVertexType(part, secondary);
256 UInt_t selectInfo = 0;
260 selectInfo = fTrackFilter->IsSelected(part);
261 if (!selectInfo) continue;
264 // add secondary tracks
265 secondary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID
273 0, // no cluster map available
276 kFALSE, // no fit performed
277 kFALSE, // no fit performed
278 AliAODTrack::kSecondary,
281 AliAODTrack* currTrack = (AliAODTrack*)tracks.Last();
282 SetChargeAndPID(part->GetPdgCode(), currTrack);
283 if (currTrack->Charge() != -99) {
284 if (currTrack->Charge() > 0) {
286 } else if (currTrack->Charge() < 0) {
291 LoopOverSecondaries(part, jTracks, jVertices, nPos, nNeg);
299 //____________________________________________________________________
300 void AliAnalysisTaskKineFilter::SetChargeAndPID(Int_t pdgCode, AliAODTrack *track) {
302 Float_t PID[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
308 PID[AliAODTrack::kUnknown] = 1.;
313 track->SetCharge(-1);
314 PID[AliAODTrack::kElectron] = 1.;
319 track->SetCharge(+1);
320 PID[AliAODTrack::kElectron] = 1.;
325 track->SetCharge(-1);
326 PID[AliAODTrack::kMuon] = 1.;
331 track->SetCharge(+1);
332 PID[AliAODTrack::kMuon] = 1.;
338 PID[AliAODTrack::kUnknown] = 1.;
343 track->SetCharge(+1);
344 PID[AliAODTrack::kPion] = 1.;
349 track->SetCharge(-1);
350 PID[AliAODTrack::kPion] = 1.;
356 PID[AliAODTrack::kUnknown] = 1.;
361 track->SetCharge(+1);
362 PID[AliAODTrack::kKaon] = 1.;
367 track->SetCharge(-1);
368 PID[AliAODTrack::kKaon] = 1.;
374 PID[AliAODTrack::kUnknown] = 1.;
379 track->SetCharge(+1);
380 PID[AliAODTrack::kProton] = 1.;
384 case -2212: // anti-p
385 track->SetCharge(-1);
386 PID[AliAODTrack::kProton] = 1.;
392 PID[AliAODTrack::kUnknown] = 1.;
398 PID[AliAODTrack::kUnknown] = 1.;
402 case -311: // anti-K0
404 PID[AliAODTrack::kUnknown] = 1.;
410 PID[AliAODTrack::kUnknown] = 1.;
416 PID[AliAODTrack::kUnknown] = 1.;
421 track->SetCharge(+1);
422 PID[AliAODTrack::kUnknown] = 1.;
427 track->SetCharge(-1);
428 PID[AliAODTrack::kUnknown] = 1.;
433 track->SetCharge(-1);
434 PID[AliAODTrack::kUnknown] = 1.;
440 PID[AliAODTrack::kUnknown] = 1.;
445 track->SetCharge(-1);
446 PID[AliAODTrack::kUnknown] = 1.;
451 track->SetCharge(-1);
452 PID[AliAODTrack::kUnknown] = 1.;
458 PID[AliAODTrack::kUnknown] = 1.;
462 case -3122: // anti-Lambda
464 PID[AliAODTrack::kUnknown] = 1.;
468 case -3222: // anti-Sigma-
469 track->SetCharge(-1);
470 PID[AliAODTrack::kUnknown] = 1.;
474 case -3212: // anti-Sigma0
476 PID[AliAODTrack::kUnknown] = 1.;
480 case -3112: // anti-Sigma+
481 track->SetCharge(+1);
482 PID[AliAODTrack::kUnknown] = 1.;
486 case -3322: // anti-Xi0
488 PID[AliAODTrack::kUnknown] = 1.;
492 case -3312: // anti-Xi+
493 track->SetCharge(+1);
496 case -3334: // anti-Omega+
497 track->SetCharge(+1);
498 PID[AliAODTrack::kUnknown] = 1.;
503 track->SetCharge(+1);
504 PID[AliAODTrack::kUnknown] = 1.;
509 track->SetCharge(-1);
510 PID[AliAODTrack::kUnknown] = 1.;
516 PID[AliAODTrack::kUnknown] = 1.;
520 case -421: // anti-D0
522 PID[AliAODTrack::kUnknown] = 1.;
527 track->SetCharge(-99);
528 PID[AliAODTrack::kUnknown] = 1.;
535 //____________________________________________________________________
536 void AliAnalysisTaskKineFilter::SetVertexType(TParticle *part, AliAODVertex *vertex)
538 // this whole thing doesn't make much sense. but anyhow...
539 AliStack* stack = MCEvent()->Stack();
540 TParticle *mother = stack->Particle(part->GetFirstMother());
541 Int_t pdgMother = mother->GetPdgCode();
542 Int_t pdgPart = part->GetPdgCode();
545 if (mother->GetNDaughters() == 2) {
546 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
547 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
549 if (!(pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
550 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 || pdgMother == 221 ||
551 TMath::Abs(pdgMother) == 3122 || TMath::Abs(pdgMother) == 3322 ||
552 pdgMother == -3212 || TMath::Abs(pdgMother) == 421 ||
553 TMath::Abs(pdgMother) == 311) // not neutral
554 && (((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
555 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
556 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
557 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
558 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // neutral
559 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
560 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
561 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
562 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
563 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) // not neutral
564 || !((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
565 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
566 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
567 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
568 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
569 && (lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
570 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
571 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
572 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
573 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)))) { // neutral
575 vertex->SetType(AliAODVertex::kKink);
581 else if (mother->GetNDaughters() == 2) {
582 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
583 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
585 if ((pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
586 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 ||
587 pdgMother == 221 || TMath::Abs(pdgMother) == 3122 ||
588 TMath::Abs(pdgMother) == 3322 || pdgMother == -3212 ||
589 TMath::Abs(pdgMother) == 421 || TMath::Abs(pdgMother) == 311) // neutral
590 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
591 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
592 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
593 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
594 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
595 && !(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
596 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
597 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
598 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
599 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) { // not neutral
601 vertex->SetType(AliAODVertex::kV0);
607 else if (mother->GetNDaughters() == 2) {
608 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
609 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
611 if ((TMath::Abs(pdgMother) == 3334 || TMath::Abs(pdgMother) == 3312 || TMath::Abs(pdgMother) == 3322) &&
612 (TMath::Abs(pdgPart) == 3122 || TMath::Abs(pdgPart) == 211 || TMath::Abs(pdgPart) == 321)
613 && ((!(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
614 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
615 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
616 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
617 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
618 && TMath::Abs(lastPdgCode) == 3122) // labmda or anti-lambda
619 || ((!(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
620 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
621 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
622 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
623 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
624 && TMath::Abs(firstPdgCode) == 3122)))) { // lambda or anti-lambda
625 vertex->SetType(AliAODVertex::kCascade);
631 else if (mother->GetNDaughters() > 2) {
633 vertex->SetType(AliAODVertex::kMulti);
638 vertex->SetType(AliAODVertex::kUndef);