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("aodTree", "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();
127 Double_t emEnergy[2] = {-999., -999.};
130 *header = AliAODHeader(aliHeader->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());
158 // Access to the AOD container of tracks
159 TClonesArray &tracks = *(aod->GetTracks());
162 aod->ResetStd(nTracks, 1);
165 for (Int_t iTrack = 0; iTrack < nPrims; ++iTrack) {
167 TParticle *part = stack->Particle(iTrack);
169 //if (part->IsPrimary()) { // this will exclude 'funny' primaries, too
170 if (kTRUE) { // this won't
171 p[0] = part->Px(); p[1] = part->Py(); p[2] = part->Pz();
172 x[0] = part->Vx(); x[1] = part->Vy(); x[2] = part->Vz();
175 // add primary vertex
176 primary = new(vertices[jVertices++])
177 AliAODVertex(x, NULL, -999., NULL, AliAODVertex::kPrimary);
180 // add primary tracks
181 primary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID,
189 0, // no ITSClusterMap
192 kFALSE, // no fit performed
193 kFALSE, // no fit preformed
194 AliAODTrack::kPrimary));
195 currTrack = (AliAODTrack*)tracks.Last();
196 SetChargeAndPID(part->GetPdgCode(), currTrack);
197 if (currTrack->Charge() != -99) {
198 if (currTrack->Charge() > 0) {
200 } else if (currTrack->Charge() < 0) {
205 LoopOverSecondaries(part);
207 } // end of track loop
209 header->SetRefMultiplicityPos(nPos);
210 header->SetRefMultiplicityNeg(nNeg);
212 cout << ":: primaries: " << nPrims << " secondaries: " << tracks.GetEntriesFast()-nPrims <<
213 " (pos: " << nPos << ", neg: " << nNeg << "), vertices: " << vertices.GetEntriesFast() <<
214 " (kinks: " << jKinks << ", V0: " << jV0s << ", cascades: " << jCascades <<
215 ", multi: " << jMultis << ")" << endl;
217 cout << " gamma: " << nGamma << " e-: " << nElectron << " e+: " << nPositron << " p: " << nProton <<
218 " pBar: " << nAntiProton << " n: " << nNeutron << " nBar: " << nAntiNeutron << " pi0: " << nPi0 <<
219 " pi+: " << nPiPlus << " pi-: " << nPiMinus << " K0: " << nK0 << " K+: " << nKPlus <<
220 " K-: " << nKMinus << endl;
222 // fill the tree for this event
224 } // end of event loop
226 aodTree->GetUserInfo()->Add(aod);
228 // write the tree to the specified file
229 outFile = aodTree->GetCurrentFile();
237 Int_t LoopOverSecondaries(TParticle *mother) {
239 if (mother->GetNDaughters() > 0) {
241 TClonesArray &vertices = *(aod->GetVertices());
242 TClonesArray &tracks = *(aod->GetTracks());
244 for (Int_t iDaughter = mother->GetFirstDaughter(); iDaughter <= mother->GetLastDaughter(); iDaughter++) {
245 TParticle *part = stack->Particle(iDaughter);
253 if (iDaughter == mother->GetFirstDaughter()) {
254 // add secondary vertex
255 secondary = new(vertices[jVertices++])
256 AliAODVertex(x, NULL, -999., tracks.Last(), AliAODVertex::kUndef);
258 SetVertexType(part, secondary);
261 // add secondary tracks
262 secondary->AddDaughter(new(tracks[jTracks++]) AliAODTrack(0, // ID
270 0, // no cluster map available
273 kFALSE, // no fit performed
274 kFALSE, // no fit performed
275 AliAODTrack::kSecondary));
277 currTrack = (AliAODTrack*)tracks.Last();
278 SetChargeAndPID(part->GetPdgCode(), currTrack);
279 if (currTrack->Charge() != -99) {
280 if (currTrack->Charge() > 0) {
282 } else if (currTrack->Charge() < 0) {
287 LoopOverSecondaries(part);
296 void SetChargeAndPID(Int_t pdgCode, AliAODTrack *track) {
298 Float_t PID[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
304 PID[AliAODTrack::kUnknown] = 1.;
310 track->SetCharge(-1);
311 PID[AliAODTrack::kElectron] = 1.;
317 track->SetCharge(+1);
318 PID[AliAODTrack::kElectron] = 1.;
324 track->SetCharge(-1);
325 PID[AliAODTrack::kMuon] = 1.;
331 track->SetCharge(+1);
332 PID[AliAODTrack::kMuon] = 1.;
339 PID[AliAODTrack::kUnknown] = 1.;
345 track->SetCharge(+1);
346 PID[AliAODTrack::kPion] = 1.;
352 track->SetCharge(-1);
353 PID[AliAODTrack::kPion] = 1.;
360 PID[AliAODTrack::kUnknown] = 1.;
366 track->SetCharge(+1);
367 PID[AliAODTrack::kKaon] = 1.;
373 track->SetCharge(-1);
374 PID[AliAODTrack::kKaon] = 1.;
381 PID[AliAODTrack::kUnknown] = 1.;
387 track->SetCharge(+1);
388 PID[AliAODTrack::kProton] = 1.;
393 case -2212: // anti-p
394 track->SetCharge(-1);
395 PID[AliAODTrack::kProton] = 1.;
402 PID[AliAODTrack::kUnknown] = 1.;
409 PID[AliAODTrack::kUnknown] = 1.;
414 case -311: // anti-K0
416 PID[AliAODTrack::kUnknown] = 1.;
423 PID[AliAODTrack::kUnknown] = 1.;
429 PID[AliAODTrack::kUnknown] = 1.;
434 track->SetCharge(+1);
435 PID[AliAODTrack::kUnknown] = 1.;
440 track->SetCharge(-1);
441 PID[AliAODTrack::kUnknown] = 1.;
446 track->SetCharge(-1);
447 PID[AliAODTrack::kUnknown] = 1.;
453 PID[AliAODTrack::kUnknown] = 1.;
458 track->SetCharge(-1);
459 PID[AliAODTrack::kUnknown] = 1.;
464 track->SetCharge(-1);
465 PID[AliAODTrack::kUnknown] = 1.;
471 PID[AliAODTrack::kUnknown] = 1.;
475 case -3122: // anti-Lambda
477 PID[AliAODTrack::kUnknown] = 1.;
481 case -3222: // anti-Sigma-
482 track->SetCharge(-1);
483 PID[AliAODTrack::kUnknown] = 1.;
487 case -3212: // anti-Sigma0
489 PID[AliAODTrack::kUnknown] = 1.;
493 case -3112: // anti-Sigma+
494 track->SetCharge(+1);
495 PID[AliAODTrack::kUnknown] = 1.;
499 case -3322: // anti-Xi0
501 PID[AliAODTrack::kUnknown] = 1.;
505 case -3312: // anti-Xi+
506 track->SetCharge(+1);
509 case -3334: // anti-Omega+
510 track->SetCharge(+1);
511 PID[AliAODTrack::kUnknown] = 1.;
516 track->SetCharge(+1);
517 PID[AliAODTrack::kUnknown] = 1.;
522 track->SetCharge(-1);
523 PID[AliAODTrack::kUnknown] = 1.;
529 PID[AliAODTrack::kUnknown] = 1.;
533 case -421: // anti-D0
535 PID[AliAODTrack::kUnknown] = 1.;
540 track->SetCharge(-99);
541 PID[AliAODTrack::kUnknown] = 1.;
549 void SetVertexType(TParticle *part, AliAODVertex *vertex) {
550 // this whole thing doesn't make much sense. but anyhow...
552 TParticle *mother = stack->Particle(part->GetFirstMother());
553 Int_t pdgMother = mother->GetPdgCode();
554 Int_t pdgPart = part->GetPdgCode();
557 if (mother->GetNDaughters() == 2) {
558 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
559 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
561 if (!(pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
562 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 || pdgMother == 221 ||
563 TMath::Abs(pdgMother) == 3122 || TMath::Abs(pdgMother) == 3322 ||
564 pdgMother == -3212 || TMath::Abs(pdgMother) == 421 ||
565 TMath::Abs(pdgMother) == 311) // not neutral
566 && (((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
567 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
568 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
569 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
570 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // neutral
571 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
572 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
573 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
574 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
575 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) // not neutral
576 || !((firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
577 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
578 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
579 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
580 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not 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)))) { // neutral
587 vertex->SetType(AliAODVertex::kKink);
593 else if (mother->GetNDaughters() == 2) {
594 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
595 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
597 if ((pdgMother == 22 || pdgMother == 111 || pdgMother == 130 ||
598 TMath::Abs(pdgMother) == 2112 || pdgMother == 310 ||
599 pdgMother == 221 || TMath::Abs(pdgMother) == 3122 ||
600 TMath::Abs(pdgMother) == 3322 || pdgMother == -3212 ||
601 TMath::Abs(pdgMother) == 421 || TMath::Abs(pdgMother) == 311) // neutral
602 && !(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
603 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
604 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
605 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
606 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
607 && !(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
608 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
609 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
610 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
611 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311)) { // not neutral
613 vertex->SetType(AliAODVertex::kV0);
619 else if (mother->GetNDaughters() == 2) {
620 Int_t firstPdgCode = stack->Particle(mother->GetFirstDaughter())->GetPdgCode();
621 Int_t lastPdgCode = stack->Particle(mother->GetLastDaughter())->GetPdgCode();
623 if ((TMath::Abs(pdgMother) == 3334 || TMath::Abs(pdgMother) == 3312 || TMath::Abs(pdgMother) == 3322) &&
624 (TMath::Abs(pdgPart) == 3122 || TMath::Abs(pdgPart) == 211 || TMath::Abs(pdgPart) == 321)
625 && ((!(firstPdgCode == 22 || firstPdgCode == 111 || firstPdgCode == 130 ||
626 TMath::Abs(firstPdgCode) == 2112 || firstPdgCode == 310 ||
627 firstPdgCode == 221 || TMath::Abs(firstPdgCode) == 3122 ||
628 TMath::Abs(firstPdgCode) == 3322 || firstPdgCode == -3212 ||
629 TMath::Abs(firstPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
630 && TMath::Abs(lastPdgCode) == 3122) // labmda or anti-lambda
631 || ((!(lastPdgCode == 22 || lastPdgCode == 111 || lastPdgCode == 130 ||
632 TMath::Abs(lastPdgCode) == 2112 || lastPdgCode == 310 ||
633 lastPdgCode == 221 || TMath::Abs(lastPdgCode) == 3122 ||
634 TMath::Abs(lastPdgCode) == 3322 || lastPdgCode == -3212 ||
635 TMath::Abs(lastPdgCode) == 421 || TMath::Abs(pdgMother) == 311) // not neutral
636 && TMath::Abs(firstPdgCode) == 3122)))) { // lambda or anti-lambda
637 vertex->SetType(AliAODVertex::kCascade);
643 else if (mother->GetNDaughters() > 2) {
645 vertex->SetType(AliAODVertex::kMulti);
650 vertex->SetType(AliAODVertex::kUndef);
653 // LocalWords: SetCharge