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;
70 //____________________________________________________________________
71 AliAnalysisTaskKineFilter& AliAnalysisTaskKineFilter::operator=(const AliAnalysisTaskKineFilter& other)
74 AliAnalysisTaskSE::operator=(other);
75 fTrackFilter = other.fTrackFilter;
79 //____________________________________________________________________
80 void AliAnalysisTaskKineFilter::UserCreateOutputObjects()
82 // Create the output container
84 OutputTree()->GetUserInfo()->Add(fTrackFilter);
88 //____________________________________________________________________
89 void AliAnalysisTaskKineFilter::Exec(Option_t */*option*/)
91 // Execute analysis for current event
94 // Fill AOD tracks from Kinematic stack
97 AliAODEvent* aod = AODEvent();
98 // aod->CreateStdContent();
100 AliStack* stack = MCEvent()->Stack();
101 Int_t nTracks = stack->GetNtrack();
102 Int_t nPrims = stack->GetNprimary();
104 AliAODVertex *primary = NULL;
112 // Access to the header
113 AliAODHeader *header = aod->GetHeader();
115 Double_t emEnergy[2] = {-999., -999.};
118 *header = AliAODHeader(MCEvent()->Header()->GetRun(),
126 -999., // muon mag. field
128 -999., // ZDCN1Energy
129 -999., // ZDCP1Energy
130 -999., // ZDCN2Energy
131 -999., // ZDCP2Energy
132 emEnergy, // emEnergy
138 // Access to the AOD container of vertices
139 TClonesArray &vertices = *(aod->GetVertices());
141 // Access to the AOD container of tracks
142 TClonesArray &tracks = *(aod->GetTracks());
144 aod->ResetStd(nTracks, 1);
147 for (Int_t iTrack = 0; iTrack < nPrims; ++iTrack) {
149 TParticle *part = stack->Particle(iTrack);
152 // add primary vertex
153 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
154 primary = new(vertices[jVertices++])
155 AliAODVertex(x, NULL, -999., NULL, AliAODVertex::kPrimary);
158 // only final particles
159 if( part->GetStatusCode() !=1 ) continue;
163 UInt_t selectInfo = 0;
165 selectInfo = fTrackFilter->IsSelected(part);
166 if (!selectInfo) continue;
169 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
170 p[0] = part->Px(); p[1] = part->Py(); p[2] = part->Pz();
172 // add primary tracks
173 primary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID,
181 0, // no ITSClusterMap
184 kFALSE, // no fit performed
185 kFALSE, // no fit preformed
186 AliAODTrack::kPrimary,
189 AliAODTrack* currTrack = (AliAODTrack*)tracks.Last();
190 SetChargeAndPID(part->GetPdgCode(), currTrack);
191 if (currTrack->Charge() != -99) {
192 if (currTrack->Charge() > 0) {
194 } else if (currTrack->Charge() < 0) {
198 LoopOverSecondaries(part, jTracks, jVertices, nPos, nNeg);
200 } // end of track loop
202 header->SetRefMultiplicityPos(nPos);
203 header->SetRefMultiplicityNeg(nNeg);
207 AliInfo(Form("primaries: %d secondaries: %d (pos: %d neg: %d), vertices: %d",
208 nPrims, tracks.GetEntriesFast()-nPrims, nPos, nNeg, vertices.GetEntriesFast() ) );
213 //____________________________________________________________________
214 Int_t AliAnalysisTaskKineFilter::LoopOverSecondaries(TParticle *mother, Int_t &jTracks, Int_t &jVertices, Int_t &nPos, Int_t &nNeg )
217 if (mother->GetNDaughters() > 0) {
219 AliStack* stack = MCEvent()->Stack();
221 TClonesArray &vertices = *(AODEvent()->GetVertices());
222 TClonesArray &tracks = *(AODEvent()->GetTracks());
225 AliAODVertex* secondary = NULL;
227 for (Int_t iDaughter = mother->GetFirstDaughter(); iDaughter <= mother->GetLastDaughter(); iDaughter++) {
228 TParticle *part = stack->Particle(iDaughter);
229 // only final particles
230 if( part->GetStatusCode() !=1 ) continue;
239 if (iDaughter == mother->GetFirstDaughter()) {
240 // add secondary vertex
241 secondary = new(vertices[jVertices++])
242 AliAODVertex(x, NULL, -999., tracks.Last(), AliAODVertex::kUndef);
244 SetVertexType(part, secondary);
247 UInt_t selectInfo = 0;
251 selectInfo = fTrackFilter->IsSelected(part);
252 if (!selectInfo) continue;
255 // add secondary tracks
256 secondary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID
264 0, // no cluster map available
267 kFALSE, // no fit performed
268 kFALSE, // no fit performed
269 AliAODTrack::kSecondary,
272 AliAODTrack* currTrack = (AliAODTrack*)tracks.Last();
273 SetChargeAndPID(part->GetPdgCode(), currTrack);
274 if (currTrack->Charge() != -99) {
275 if (currTrack->Charge() > 0) {
277 } else if (currTrack->Charge() < 0) {
282 LoopOverSecondaries(part, jTracks, jVertices, nPos, nNeg);
290 //____________________________________________________________________
291 void AliAnalysisTaskKineFilter::SetChargeAndPID(Int_t pdgCode, AliAODTrack *track) {
293 Float_t PID[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
299 PID[AliAODTrack::kUnknown] = 1.;
304 track->SetCharge(-1);
305 PID[AliAODTrack::kElectron] = 1.;
310 track->SetCharge(+1);
311 PID[AliAODTrack::kElectron] = 1.;
316 track->SetCharge(-1);
317 PID[AliAODTrack::kMuon] = 1.;
322 track->SetCharge(+1);
323 PID[AliAODTrack::kMuon] = 1.;
329 PID[AliAODTrack::kUnknown] = 1.;
334 track->SetCharge(+1);
335 PID[AliAODTrack::kPion] = 1.;
340 track->SetCharge(-1);
341 PID[AliAODTrack::kPion] = 1.;
347 PID[AliAODTrack::kUnknown] = 1.;
352 track->SetCharge(+1);
353 PID[AliAODTrack::kKaon] = 1.;
358 track->SetCharge(-1);
359 PID[AliAODTrack::kKaon] = 1.;
365 PID[AliAODTrack::kUnknown] = 1.;
370 track->SetCharge(+1);
371 PID[AliAODTrack::kProton] = 1.;
375 case -2212: // anti-p
376 track->SetCharge(-1);
377 PID[AliAODTrack::kProton] = 1.;
383 PID[AliAODTrack::kUnknown] = 1.;
389 PID[AliAODTrack::kUnknown] = 1.;
393 case -311: // anti-K0
395 PID[AliAODTrack::kUnknown] = 1.;
401 PID[AliAODTrack::kUnknown] = 1.;
407 PID[AliAODTrack::kUnknown] = 1.;
412 track->SetCharge(+1);
413 PID[AliAODTrack::kUnknown] = 1.;
418 track->SetCharge(-1);
419 PID[AliAODTrack::kUnknown] = 1.;
424 track->SetCharge(-1);
425 PID[AliAODTrack::kUnknown] = 1.;
431 PID[AliAODTrack::kUnknown] = 1.;
436 track->SetCharge(-1);
437 PID[AliAODTrack::kUnknown] = 1.;
442 track->SetCharge(-1);
443 PID[AliAODTrack::kUnknown] = 1.;
449 PID[AliAODTrack::kUnknown] = 1.;
453 case -3122: // anti-Lambda
455 PID[AliAODTrack::kUnknown] = 1.;
459 case -3222: // anti-Sigma-
460 track->SetCharge(-1);
461 PID[AliAODTrack::kUnknown] = 1.;
465 case -3212: // anti-Sigma0
467 PID[AliAODTrack::kUnknown] = 1.;
471 case -3112: // anti-Sigma+
472 track->SetCharge(+1);
473 PID[AliAODTrack::kUnknown] = 1.;
477 case -3322: // anti-Xi0
479 PID[AliAODTrack::kUnknown] = 1.;
483 case -3312: // anti-Xi+
484 track->SetCharge(+1);
487 case -3334: // anti-Omega+
488 track->SetCharge(+1);
489 PID[AliAODTrack::kUnknown] = 1.;
494 track->SetCharge(+1);
495 PID[AliAODTrack::kUnknown] = 1.;
500 track->SetCharge(-1);
501 PID[AliAODTrack::kUnknown] = 1.;
507 PID[AliAODTrack::kUnknown] = 1.;
511 case -421: // anti-D0
513 PID[AliAODTrack::kUnknown] = 1.;
518 track->SetCharge(-99);
519 PID[AliAODTrack::kUnknown] = 1.;
526 //____________________________________________________________________
527 void AliAnalysisTaskKineFilter::SetVertexType(TParticle *part, AliAODVertex *vertex)
529 // this whole thing doesn't make much sense. but anyhow...
530 AliStack* stack = MCEvent()->Stack();
531 TParticle *mother = stack->Particle(part->GetFirstMother());
532 Int_t pdgMother = mother->GetPdgCode();
533 Int_t pdgPart = part->GetPdgCode();
536 if (mother->GetNDaughters() == 2) {
537 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
538 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
540 if (!(pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
541 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 || pdgMother == 221 ||
542 TMath::Abs(pdgMother) == 3122 || TMath::Abs(pdgMother) == 3322 ||
543 pdgMother == -3212 || TMath::Abs(pdgMother) == 421 ||
544 TMath::Abs(pdgMother) == 311) // not neutral
545 && (((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
546 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
547 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
548 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
549 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // neutral
550 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
551 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
552 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
553 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
554 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) // not neutral
555 || !((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
556 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
557 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
558 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
559 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
560 && (lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
561 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
562 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
563 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
564 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)))) { // neutral
566 vertex->SetType(AliAODVertex::kKink);
572 else if (mother->GetNDaughters() == 2) {
573 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
574 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
576 if ((pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
577 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 ||
578 pdgMother == 221 || TMath::Abs(pdgMother) == 3122 ||
579 TMath::Abs(pdgMother) == 3322 || pdgMother == -3212 ||
580 TMath::Abs(pdgMother) == 421 || TMath::Abs(pdgMother) == 311) // neutral
581 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
582 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
583 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
584 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
585 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
586 && !(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
587 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
588 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
589 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
590 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) { // not neutral
592 vertex->SetType(AliAODVertex::kV0);
598 else if (mother->GetNDaughters() == 2) {
599 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
600 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
602 if ((TMath::Abs(pdgMother) == 3334 || TMath::Abs(pdgMother) == 3312 || TMath::Abs(pdgMother) == 3322) &&
603 (TMath::Abs(pdgPart) == 3122 || TMath::Abs(pdgPart) == 211 || TMath::Abs(pdgPart) == 321)
604 && ((!(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
605 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
606 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
607 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
608 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
609 && TMath::Abs(lastPdgCode) == 3122) // labmda or anti-lambda
610 || ((!(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
611 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
612 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
613 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
614 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
615 && TMath::Abs(firstPdgCode) == 3122)))) { // lambda or anti-lambda
616 vertex->SetType(AliAODVertex::kCascade);
622 else if (mother->GetNDaughters() > 2) {
624 vertex->SetType(AliAODVertex::kMulti);
629 vertex->SetType(AliAODVertex::kUndef);