1 #if !defined(__CINT__) || defined(__MAKECINT__)
9 #include "AliRunLoader.h"
10 #include "AliHeader.h"
12 #include "TParticle.h"
14 #include "AliAODEvent.h"
15 #include "AliAODHeader.h"
16 #include "AliAODVertex.h"
17 #include "AliAODTrack.h"
18 #include "AliAODCluster.h"
22 // function declarations
23 void SetVertexType(TParticle *part, AliAODVertex *vertex);
24 void SetChargeAndPID(Int_t pdgCode, AliAODTrack *track);
25 Int_t LoopOverSecondaries(TParticle *mother);
26 void CreateAODfromKineTree(const char *inFileName = "galice.root",
27 const char *outFileName = "AliAOD.root");
30 AliAODEvent *aod = NULL;
31 AliStack *stack = NULL;
48 Int_t nAntiProton = 0;
50 Int_t nAntiNeutron = 0;
58 // global arrays and pointers
61 Float_t *cov = NULL; // set to NULL because not provided
62 Float_t *pid = NULL; // set to NULL because not provided
63 AliAODVertex *primary = NULL;
64 AliAODVertex *secondary = NULL;
65 AliAODTrack *currTrack = NULL;
67 void CreateAODfromKineTree(const char *inFileName,
68 const char *outFileName) {
69 printf("This macro works only correctly in comiled mode!\n");
70 /* This probem is due to CINT which gets confused if arrays are reused */
72 // create an AliAOD object
73 aod = new AliAODEvent();
74 aod->CreateStdContent();
77 TFile *outFile = TFile::Open(outFileName, "RECREATE");
80 TTree *aodTree = new TTree("AOD", "AliAOD tree");
81 aodTree->Branch(aod->GetList());
83 AliRunLoader *runLoader;
86 runLoader = AliRunLoader::Open(inFileName);
87 if (!runLoader) return;
89 runLoader->LoadHeader();
90 runLoader->LoadKinematics();
92 Int_t nEvents = runLoader->GetNumberOfEvents();
94 // loop over events and fill them
95 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
97 cout << "Event " << iEvent+1 << "/" << nEvents;
98 runLoader->GetEvent(iEvent);
100 AliHeader *aliHeader = runLoader->GetHeader();
101 stack = runLoader->Stack();
102 Int_t nTracks = stack->GetNtrack();
103 Int_t nPrims = stack->GetNprimary();
124 // Access to the header
125 AliAODHeader *header = aod->GetHeader();
128 *header = AliAODHeader(aliHeader->GetRun(),
136 -999., // muon mag. field
147 // Access to the AOD container of vertices
148 TClonesArray &vertices = *(aod->GetVertices());
155 // Access to the AOD container of tracks
156 TClonesArray &tracks = *(aod->GetTracks());
159 aod->ResetStd(nTracks, 1);
162 for (Int_t iTrack = 0; iTrack < nPrims; ++iTrack) {
164 TParticle *part = stack->Particle(iTrack);
166 //if (part->IsPrimary()) { // this will exclude 'funny' primaries, too
167 if (kTRUE) { // this won't
168 p[0] = part->Px(); p[1] = part->Py(); p[2] = part->Pz();
169 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
172 // add primary vertex
173 primary = new(vertices[jVertices++])
174 AliAODVertex(x, NULL, -999., NULL, AliAODVertex::kPrimary);
177 // add primary tracks
178 primary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID,
186 0, // no ITSClusterMap
189 kFALSE, // no fit performed
190 kFALSE, // no fit preformed
191 AliAODTrack::kPrimary));
192 currTrack = (AliAODTrack*)tracks.Last();
193 SetChargeAndPID(part->GetPdgCode(), currTrack);
194 if (currTrack->Charge() != -99) {
195 if (currTrack->Charge() > 0) {
197 } else if (currTrack->Charge() < 0) {
202 LoopOverSecondaries(part);
204 } // end of track loop
206 header->SetRefMultiplicityPos(nPos);
207 header->SetRefMultiplicityNeg(nNeg);
209 cout << ":: primaries: " << nPrims << " secondaries: " << tracks.GetEntriesFast()-nPrims <<
210 " (pos: " << nPos << ", neg: " << nNeg << "), vertices: " << vertices.GetEntriesFast() <<
211 " (kinks: " << jKinks << ", V0: " << jV0s << ", cascades: " << jCascades <<
212 ", multi: " << jMultis << ")" << endl;
214 cout << " gamma: " << nGamma << " e-: " << nElectron << " e+: " << nPositron << " p: " << nProton <<
215 " pBar: " << nAntiProton << " n: " << nNeutron << " nBar: " << nAntiNeutron << " pi0: " << nPi0 <<
216 " pi+: " << nPiPlus << " pi-: " << nPiMinus << " K0: " << nK0 << " K+: " << nKPlus <<
217 " K-: " << nKMinus << endl;
219 // fill the tree for this event
221 } // end of event loop
223 aodTree->GetUserInfo()->Add(aod);
225 // write the tree to the specified file
226 outFile = aodTree->GetCurrentFile();
234 Int_t LoopOverSecondaries(TParticle *mother) {
236 if (mother->GetNDaughters() > 0) {
238 TClonesArray &vertices = *(aod->GetVertices());
239 TClonesArray &tracks = *(aod->GetTracks());
241 for (Int_t iDaughter = mother->GetFirstDaughter(); iDaughter <= mother->GetLastDaughter(); iDaughter++) {
242 TParticle *part = stack->Particle(iDaughter);
250 if (iDaughter == mother->GetFirstDaughter()) {
251 // add secondary vertex
252 secondary = new(vertices[jVertices++])
253 AliAODVertex(x, NULL, -999., tracks.Last(), AliAODVertex::kUndef);
255 SetVertexType(part, secondary);
258 // add secondary tracks
259 secondary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID
267 0, // no cluster map available
270 kFALSE, // no fit performed
271 kFALSE, // no fit performed
272 AliAODTrack::kSecondary));
274 currTrack = (AliAODTrack*)tracks.Last();
275 SetChargeAndPID(part->GetPdgCode(), currTrack);
276 if (currTrack->Charge() != -99) {
277 if (currTrack->Charge() > 0) {
279 } else if (currTrack->Charge() < 0) {
284 LoopOverSecondaries(part);
293 void SetChargeAndPID(Int_t pdgCode, AliAODTrack *track) {
295 Float_t PID[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
301 PID[AliAODTrack::kUnknown] = 1.;
307 track->SetCharge(-1);
308 PID[AliAODTrack::kElectron] = 1.;
314 track->SetCharge(+1);
315 PID[AliAODTrack::kElectron] = 1.;
321 track->SetCharge(-1);
322 PID[AliAODTrack::kMuon] = 1.;
328 track->SetCharge(+1);
329 PID[AliAODTrack::kMuon] = 1.;
336 PID[AliAODTrack::kUnknown] = 1.;
342 track->SetCharge(+1);
343 PID[AliAODTrack::kPion] = 1.;
349 track->SetCharge(-1);
350 PID[AliAODTrack::kPion] = 1.;
357 PID[AliAODTrack::kUnknown] = 1.;
363 track->SetCharge(+1);
364 PID[AliAODTrack::kKaon] = 1.;
370 track->SetCharge(-1);
371 PID[AliAODTrack::kKaon] = 1.;
378 PID[AliAODTrack::kUnknown] = 1.;
384 track->SetCharge(+1);
385 PID[AliAODTrack::kProton] = 1.;
390 case -2212: // anti-p
391 track->SetCharge(-1);
392 PID[AliAODTrack::kProton] = 1.;
399 PID[AliAODTrack::kUnknown] = 1.;
406 PID[AliAODTrack::kUnknown] = 1.;
411 case -311: // anti-K0
413 PID[AliAODTrack::kUnknown] = 1.;
420 PID[AliAODTrack::kUnknown] = 1.;
426 PID[AliAODTrack::kUnknown] = 1.;
431 track->SetCharge(+1);
432 PID[AliAODTrack::kUnknown] = 1.;
437 track->SetCharge(-1);
438 PID[AliAODTrack::kUnknown] = 1.;
443 track->SetCharge(-1);
444 PID[AliAODTrack::kUnknown] = 1.;
450 PID[AliAODTrack::kUnknown] = 1.;
455 track->SetCharge(-1);
456 PID[AliAODTrack::kUnknown] = 1.;
461 track->SetCharge(-1);
462 PID[AliAODTrack::kUnknown] = 1.;
468 PID[AliAODTrack::kUnknown] = 1.;
472 case -3122: // anti-Lambda
474 PID[AliAODTrack::kUnknown] = 1.;
478 case -3222: // anti-Sigma-
479 track->SetCharge(-1);
480 PID[AliAODTrack::kUnknown] = 1.;
484 case -3212: // anti-Sigma0
486 PID[AliAODTrack::kUnknown] = 1.;
490 case -3112: // anti-Sigma+
491 track->SetCharge(+1);
492 PID[AliAODTrack::kUnknown] = 1.;
496 case -3322: // anti-Xi0
498 PID[AliAODTrack::kUnknown] = 1.;
502 case -3312: // anti-Xi+
503 track->SetCharge(+1);
506 case -3334: // anti-Omega+
507 track->SetCharge(+1);
508 PID[AliAODTrack::kUnknown] = 1.;
513 track->SetCharge(+1);
514 PID[AliAODTrack::kUnknown] = 1.;
519 track->SetCharge(-1);
520 PID[AliAODTrack::kUnknown] = 1.;
526 PID[AliAODTrack::kUnknown] = 1.;
530 case -421: // anti-D0
532 PID[AliAODTrack::kUnknown] = 1.;
537 track->SetCharge(-99);
538 PID[AliAODTrack::kUnknown] = 1.;
546 void SetVertexType(TParticle *part, AliAODVertex *vertex) {
547 // this whole thing doesn't make much sense. but anyhow...
549 TParticle *mother = stack->Particle(part->GetFirstMother());
550 Int_t pdgMother = mother->GetPdgCode();
551 Int_t pdgPart = part->GetPdgCode();
554 if (mother->GetNDaughters() == 2) {
555 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
556 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
558 if (!(pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
559 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 || pdgMother == 221 ||
560 TMath::Abs(pdgMother) == 3122 || TMath::Abs(pdgMother) == 3322 ||
561 pdgMother == -3212 || TMath::Abs(pdgMother) == 421 ||
562 TMath::Abs(pdgMother) == 311) // not neutral
563 && (((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
564 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
565 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
566 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
567 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // neutral
568 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
569 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
570 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
571 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
572 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) // not neutral
573 || !((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
574 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
575 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
576 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
577 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
578 && (lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
579 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
580 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
581 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
582 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)))) { // neutral
584 vertex->SetType(AliAODVertex::kKink);
590 else if (mother->GetNDaughters() == 2) {
591 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
592 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
594 if ((pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
595 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 ||
596 pdgMother == 221 || TMath::Abs(pdgMother) == 3122 ||
597 TMath::Abs(pdgMother) == 3322 || pdgMother == -3212 ||
598 TMath::Abs(pdgMother) == 421 || TMath::Abs(pdgMother) == 311) // neutral
599 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
600 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
601 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
602 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
603 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
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
610 vertex->SetType(AliAODVertex::kV0);
616 else if (mother->GetNDaughters() == 2) {
617 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
618 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
620 if ((TMath::Abs(pdgMother) == 3334 || TMath::Abs(pdgMother) == 3312 || TMath::Abs(pdgMother) == 3322) &&
621 (TMath::Abs(pdgPart) == 3122 || TMath::Abs(pdgPart) == 211 || TMath::Abs(pdgPart) == 321)
622 && ((!(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
623 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
624 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
625 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
626 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
627 && TMath::Abs(lastPdgCode) == 3122) // labmda or anti-lambda
628 || ((!(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
629 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
630 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
631 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
632 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
633 && TMath::Abs(firstPdgCode) == 3122)))) { // lambda or anti-lambda
634 vertex->SetType(AliAODVertex::kCascade);
640 else if (mother->GetNDaughters() > 2) {
642 vertex->SetType(AliAODVertex::kMulti);
647 vertex->SetType(AliAODVertex::kUndef);
650 // LocalWords: SetCharge