New ITS geometry class which uses direct queries to TGeoManager
[u/mrichter/AliRoot.git] / STEER / CreateAODfromESD.C
CommitLineData
cf22b3fc 1#include <Riostream.h>
2#include <TFile.h>
3#include <TTree.h>
4#include <TMath.h>
df9db588 5
6#include "AliAODEvent.h"
cf22b3fc 7#include "AliAODVertex.h"
8#include "AliAODTrack.h"
9#include "AliAODNeutral.h"
10
df9db588 11#include "AliESD.h"
12#include "AliESDtrack.h"
13#include "AliESDVertex.h"
14#include "AliESDv0.h"
cf22b3fc 15#include "AliESDCascade.h"
16#include "AliESDCaloCluster.h"
df9db588 17
cf22b3fc 18void CreateAODfromESD(const char *inFileName = "AliESDs.root",
19 const char *outFileName = "AliAOD.root") {
df9db588 20
21 // create an AliAOD object
22 AliAODEvent *aod = new AliAODEvent();
23 aod->CreateStdContent();
24
25 // open the file
26 TFile *outFile = TFile::Open(outFileName, "RECREATE");
27
28 // create the tree
29 TTree *aodTree = new TTree("AOD", "AliAOD tree");
30 aodTree->Branch(aod->GetList());
31
32 // connect to ESD
33 TFile *inFile = TFile::Open(inFileName, "READ");
34 TTree *t = (TTree*) inFile->Get("esdTree");
35 TBranch *b = t->GetBranch("ESD");
36 AliESD *esd = 0;
37 b->SetAddress(&esd);
38
39 Int_t nEvents = b->GetEntries();
40
41 // loop over events and fill them
42 for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
43 b->GetEntry(iEvent);
44
45 // Multiplicity information needed by the header (to be revised!)
46 Int_t nTracks = esd->GetNumberOfTracks();
47 Int_t nPosTracks = 0;
48 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack)
49 if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
50
51 // create the header
52 aod->AddHeader(new AliAODHeader(esd ->GetEventNumber(),
53 esd->GetRunNumber(),
54 nTracks,
55 nPosTracks,
56 nTracks-nPosTracks,
57 esd->GetMagneticField(),
58 -999., // centrality; to be filled, still
59 esd->GetTriggerMask(),
60 esd->GetTriggerCluster(),
61 esd->GetEventType()));
62
63 Int_t nV0s = esd->GetNumberOfV0s();
64 Int_t nCascades = esd->GetNumberOfCascades();
65 Int_t nKinks = esd->GetNumberOfKinks();
66 Int_t nVertices = nV0s + nCascades + nKinks;
67
68 aod->ResetStd(nTracks, nVertices);
69
70
71 // Array to take into account the tracks already added to the AOD
72 Bool_t * usedTrack = NULL;
73 if (nTracks>0) {
74 usedTrack = new Bool_t[nTracks];
75 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
76 }
77 // Array to take into account the V0s already added to the AOD
78 Bool_t * usedV0 = NULL;
79 if (nV0s>0) {
80 usedV0 = new Bool_t[nV0s];
81 for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
82 }
cf22b3fc 83 // Array to take into account the kinks already added to the AOD
84 Bool_t * usedKink = NULL;
85 if (nKinks>0) {
86 usedKink = new Bool_t[nKinks];
87 for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
88 }
df9db588 89
90 // Access to the AOD container of vertices
91 TClonesArray &vertices = *(aod->GetVertices());
92 Int_t jVertices=0;
93
94 // Access to the AOD container of tracks
95 TClonesArray &tracks = *(aod->GetTracks());
96 Int_t jTracks=0;
97
98 // Add primary vertex. The primary tracks will be defined
99 // after the loops on the composite objects (V0, cascades, kinks)
100 const AliESDVertex *vtx = esd->GetPrimaryVertex();
101
102 Double_t pos[3];
103 vtx->GetXYZ(pos); // position
104 Double_t cov[6];
105 vtx->GetCovMatrix(cov); //covariance matrix
106
107 AliAODVertex * primary = new(vertices[jVertices++])
108 AliAODVertex(pos, cov, vtx->GetChi2(), NULL, AliAODVertex::kPrimary);
109
110 // Create vertices starting from the most complex objects
111
112 // Cascades
113 for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
114 AliESDcascade *cascade = esd->GetCascade(nCascade);
115
116 Double_t posXi[3];
117 cascade->GetXYZ(posXi[0], posXi[1], posXi[2]);
118 Double_t covXi[6];
119 cascade->GetPosCovXi(covXi);
120
121 // Add the cascade vertex
122 AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(posXi,
123 covXi,
124 cascade->GetChi2Xi(),
125 primary,
126 AliAODVertex::kCascade);
127
128 primary->AddDaughter(vcascade);
129
130 // Add the V0 from the cascade. The ESD class have to be optimized...
131 // Now we have to search for the corresponding Vo in the list of V0s
132 // using the indeces of the positive and negative tracks
133
134 Int_t posFromV0 = cascade->GetPindex();
135 Int_t negFromV0 = cascade->GetNindex();
136
137
138 AliESDv0 * v0 = 0x0;
139 Int_t indV0 = -1;
140
141 for (Int_t iV0=0; iV0<nV0s; ++iV0) {
142
143 v0 = esd->GetV0(iV0);
144 Int_t pos = v0->GetPindex();
145 Int_t neg = v0->GetNindex();
146
147 if (pos==posFromV0 && neg==negFromV0) {
148 indV0 = iV0;
149 break;
150 }
151 }
152
153 AliAODVertex * vV0FromCascade = 0x0;
154
cf22b3fc 155 if (indV0>-1 && !usedV0[indV0] ) {
df9db588 156
cf22b3fc 157 // the V0 exists in the array of V0s and is not used
158
df9db588 159 usedV0[indV0] = kTRUE;
160
161 Double_t posV0[3];
162 v0->GetXYZ(posV0[0], posV0[1], posV0[2]);
163 Double_t covV0[6];
164 v0->GetPosCov(covV0);
cf22b3fc 165
df9db588 166 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(posV0,
167 covV0,
168 v0->GetChi2V0(),
169 vcascade,
170 AliAODVertex::kV0);
171 } else {
172
cf22b3fc 173 // the V0 doesn't exist in the array of V0s or was used
174 cerr << "Error: event " << iEvent << " cascade " << nCascade
175 << " The V0 " << indV0
176 << " doesn't exist in the array of V0s or was used!" << endl;
177
df9db588 178 Double_t posV0[3];
179 cascade->GetXYZ(posV0[0], posV0[1], posV0[2]);
180 Double_t covV0[6];
181 cascade->GetPosCov(covV0);
182
183 vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(posV0,
184 covV0,
185 v0->GetChi2V0(),
186 vcascade,
187 AliAODVertex::kV0);
188 vcascade->AddDaughter(vV0FromCascade);
189 }
190
191 // Add the positive tracks from the V0
192
cf22b3fc 193 if (! usedTrack[posFromV0]) {
194
195 usedTrack[posFromV0] = kTRUE;
df9db588 196
197 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
cf22b3fc 198
df9db588 199 Double_t p[3];
200 esdTrack->GetPxPyPz(p);
201
202 Double_t x[3];
203 esdTrack->GetXYZ(x);
204
205 Double_t cov[21];
206 esdTrack->GetCovarianceXYZPxPyPz(cov);
cf22b3fc 207
df9db588 208 Double_t pid[10];
209 esdTrack->GetESDpid(pid);
210
211 vV0FromCascade->AddDaughter(
cf22b3fc 212 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
df9db588 213 esdTrack->GetLabel(),
214 p, kTRUE,
215 x,
216 kFALSE,
217 cov,
218 (Short_t)esdTrack->GetSign(),
219 esdTrack->GetITSClusterMap(),
220 pid,
221 vV0FromCascade,
222 AliAODTrack::kSecondary)
223 );
224 }
cf22b3fc 225 else {
226 cerr << "Error: event " << iEvent << " cascade " << nCascade
227 << " track " << posFromV0 << " has already been used!" << endl;
228 }
df9db588 229
230 // Add the negative tracks from the V0
231
cf22b3fc 232 if (!usedTrack[negFromV0]) {
233
234 usedTrack[negFromV0] = kTRUE;
235
df9db588 236 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
cf22b3fc 237
df9db588 238 Double_t p[3];
239 esdTrack->GetPxPyPz(p);
240
241 Double_t x[3];
242 esdTrack->GetXYZ(x);
243
244 Double_t cov[21];
245 esdTrack->GetCovarianceXYZPxPyPz(cov);
246
247 Double_t pid[10];
248 esdTrack->GetESDpid(pid);
cf22b3fc 249
df9db588 250 vV0FromCascade->AddDaughter(
251 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
252 esdTrack->GetLabel(),
253 p,
254 kTRUE,
255 x,
256 kFALSE,
257 cov,
258 (Short_t)esdTrack->GetSign(),
259 esdTrack->GetITSClusterMap(),
260 pid,
261 vV0FromCascade,
262 AliAODTrack::kSecondary)
263 );
264 }
cf22b3fc 265 else {
266 cerr << "Error: event " << iEvent << " cascade " << nCascade
267 << " track " << negFromV0 << " has already been used!" << endl;
268 }
df9db588 269
270 // Add the bachelor track from the cascade
271
272 Int_t bachelor = cascade->GetBindex();
df9db588 273
cf22b3fc 274 if(!usedTrack[bachelor]) {
275
276 usedTrack[bachelor] = kTRUE;
277
278 AliESDtrack *esdTrack = esd->GetTrack(bachelor);
279
df9db588 280 Double_t p[3];
281 esdTrack->GetPxPyPz(p);
282
283 Double_t x[3];
284 esdTrack->GetXYZ(x);
285
286 Double_t cov[21];
287 esdTrack->GetCovarianceXYZPxPyPz(cov);
288
289 Double_t pid[10];
290 esdTrack->GetESDpid(pid);
291
292 vcascade->AddDaughter(
293 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
294 esdTrack->GetLabel(),
295 p,
296 kTRUE,
297 x,
298 kFALSE,
299 cov,
300 (Short_t)esdTrack->GetSign(),
301 esdTrack->GetITSClusterMap(),
302 pid,
303 vcascade,
304 AliAODTrack::kSecondary)
305 );
306 }
cf22b3fc 307 else {
308 cerr << "Error: event " << iEvent << " cascade " << nCascade
309 << " track " << bachelor << " has already been used!" << endl;
310 }
df9db588 311
312 // Add the primary track of the cascade (if any)
313
cf22b3fc 314 } // end of the loop on cascades
df9db588 315
316 // V0s
317
318 for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
319
320 if (usedV0[nV0]) continue; // skip if aready added to the AOD
321
322 AliESDv0 *v0 = esd->GetV0(nV0);
323
324 Double_t posV0[3];
325 v0->GetXYZ(posV0[0], posV0[1], posV0[2]);
326 Double_t covV0[6];
327 v0->GetPosCov(covV0);
328
329 AliAODVertex * vV0 =
330 new(vertices[jVertices++]) AliAODVertex(posV0,
331 covV0,
332 v0->GetChi2V0(),
333 primary,
334 AliAODVertex::kV0);
335 primary->AddDaughter(vV0);
336
337 Int_t posFromV0 = v0->GetPindex();
338 Int_t negFromV0 = v0->GetNindex();
339
340 // Add the positive tracks from the V0
341
cf22b3fc 342 if (!usedTrack[posFromV0]) {
343
344 usedTrack[posFromV0] = kTRUE;
df9db588 345
346 AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
347
348 Double_t p[3];
349 esdTrack->GetPxPyPz(p);
350
351 Double_t x[3];
352 esdTrack->GetXYZ(x);
353
354 Double_t cov[21];
355 esdTrack->GetCovarianceXYZPxPyPz(cov);
cf22b3fc 356
df9db588 357 Double_t pid[10];
358 esdTrack->GetESDpid(pid);
359
360 vV0->AddDaughter(
361 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
362 esdTrack->GetLabel(),
363 p, kTRUE,
364 x,
365 kFALSE,
366 cov,
367 (Short_t)esdTrack->GetSign(),
368 esdTrack->GetITSClusterMap(),
369 pid,
370 vV0,
371 AliAODTrack::kSecondary)
372 );
373 }
cf22b3fc 374 else {
375 cerr << "Error: event " << iEvent << " V0 " << nV0
376 << " track " << posFromV0 << " has already been used!" << endl;
377 }
df9db588 378
379 // Add the negative tracks from the V0
380
cf22b3fc 381 if (!usedTrack[negFromV0]) {
382
383 usedTrack[negFromV0] = kTRUE;
df9db588 384
385 AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
386
387 Double_t p[3];
388 esdTrack->GetPxPyPz(p);
389
390 Double_t x[3];
391 esdTrack->GetXYZ(x);
392
393 Double_t cov[21];
394 esdTrack->GetCovarianceXYZPxPyPz(cov);
395
396 Double_t pid[10];
397 esdTrack->GetESDpid(pid);
398
399 vV0->AddDaughter(
400 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
401 esdTrack->GetLabel(),
402 p,
403 kTRUE,
404 x,
405 kFALSE,
406 cov,
407 (Short_t)esdTrack->GetSign(),
408 esdTrack->GetITSClusterMap(),
409 pid,
410 vV0,
411 AliAODTrack::kSecondary)
412 );
413 }
cf22b3fc 414 else {
415 cerr << "Error: event " << iEvent << " V0 " << nV0
416 << " track " << negFromV0 << " has already been used!" << endl;
417 }
df9db588 418
cf22b3fc 419 } // end of the loop on V0s
df9db588 420
cf22b3fc 421 // Kinks: it is a big mess the access to the information in the kinks
422 // The loop is on the tracks in order to find the mother and daugther of each kink
423
424
425 for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
426
427
428 AliESDtrack * track = esd->GetTrack(iTrack);
429
430
431 Int_t ikink = track->GetKinkIndex(0);
432
433 if (ikink) {
434 // Negative kink index: mother, positive: daughter
435
436 // Search for the second track of the kink
437
438 for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
439
440 AliESDtrack * track1 = esd->GetTrack(jTrack);
441
442 Int_t jkink = track1->GetKinkIndex(0);
443
444 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
445
446 // The two tracks are from the same kink
447
448 if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
449
450 Int_t imother = -1;
451 Int_t idaughter = -1;
452
453 if (ikink<0 && jkink>0) {
454
455 imother = iTrack;
456 idaughter = jTrack;
457 }
458 else if (ikink>0 && jkink<0) {
459
460 imother = jTrack;
461 idaughter = iTrack;
462 }
463 else {
464 cerr << "Error: Wrong combination of kink indexes: "
465 << ikink << " " << jkink << endl;
466 continue;
467 }
468
469 // Add the mother track
470
471 AliAODTrack * mother = NULL;
472
473 if (!usedTrack[imother]) {
474
475 usedTrack[imother] = kTRUE;
476
477 AliESDtrack *esdTrack = esd->GetTrack(imother);
478
479 Double_t p[3];
480 esdTrack->GetPxPyPz(p);
481
482 Double_t x[3];
483 esdTrack->GetXYZ(x);
484
485 Double_t cov[21];
486 esdTrack->GetCovarianceXYZPxPyPz(cov);
487
488 Double_t pid[10];
489 esdTrack->GetESDpid(pid);
490
491
492 mother =
493 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
494 esdTrack->GetLabel(),
495 p,
496 kTRUE,
497 x,
498 kFALSE,
499 cov,
500 (Short_t)esdTrack->GetSign(),
501 esdTrack->GetITSClusterMap(),
502 pid,
503 primary,
504 AliAODTrack::kPrimary);
505 primary->AddDaughter(mother);
506 }
507 else {
508 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
509 << " track " << imother << " has already been used!" << endl;
510 }
511
512 // Add the kink vertex
513 AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
514
515 AliAODVertex * vkink =
516 new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
517 NULL,
518 0.,
519 mother,
520 AliAODVertex::kKink);
521 // Add the daughter track
522
523 AliAODTrack * daughter = NULL;
524
525 if (!usedTrack[idaughter]) {
526
527 usedTrack[idaughter] = kTRUE;
528
529 AliESDtrack *esdTrack = esd->GetTrack(idaughter);
530
531 Double_t p[3];
532 esdTrack->GetPxPyPz(p);
533
534 Double_t x[3];
535 esdTrack->GetXYZ(x);
536
537 Double_t cov[21];
538 esdTrack->GetCovarianceXYZPxPyPz(cov);
539
540 Double_t pid[10];
541 esdTrack->GetESDpid(pid);
542
543
544 daughter =
545 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
546 esdTrack->GetLabel(),
547 p,
548 kTRUE,
549 x,
550 kFALSE,
551 cov,
552 (Short_t)esdTrack->GetSign(),
553 esdTrack->GetITSClusterMap(),
554 pid,
555 vkink,
556 AliAODTrack::kPrimary);
557 vkink->AddDaughter(daughter);
558 }
559 else {
560 cerr << "Error: event " << iEvent << " kink " << TMath::Abs(ikink)-1
561 << " track " << idaughter << " has already been used!" << endl;
562 }
563
564
565 }
566 }
567
568 }
df9db588 569
df9db588 570 }
cf22b3fc 571
df9db588 572
cf22b3fc 573 // Tracks (primary and orphan)
df9db588 574
575 for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
576
577
578 if (usedTrack[nTrack]) continue;
579
580 AliESDtrack *esdTrack = esd->GetTrack(nTrack);
581
582 Double_t p[3];
583 esdTrack->GetPxPyPz(p);
584
585 Double_t x[3];
586 esdTrack->GetXYZ(x);
587
588 Double_t cov[21];
589 esdTrack->GetCovarianceXYZPxPyPz(cov);
590
591 Double_t pid[10];
592 esdTrack->GetESDpid(pid);
cf22b3fc 593
594 Float_t impactXY, impactZ;
595
596 esdTrack->GetImpactParameters(impactXY,impactZ);
597
598 if (impactXY<3) {
599 // track inside the beam pipe
df9db588 600
cf22b3fc 601 primary->AddDaughter(
602 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
df9db588 603 esdTrack->GetLabel(),
604 p,
605 kTRUE,
606 x,
607 kFALSE,
608 cov,
609 (Short_t)esdTrack->GetSign(),
610 esdTrack->GetITSClusterMap(),
611 pid,
612 primary,
613 AliAODTrack::kPrimary)
cf22b3fc 614 );
615 }
616 else {
617 // outside the beam pipe: orphan track
618 new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
619 esdTrack->GetLabel(),
620 p,
621 kTRUE,
622 x,
623 kFALSE,
624 cov,
625 (Short_t)esdTrack->GetSign(),
626 esdTrack->GetITSClusterMap(),
627 pid,
628 NULL,
629 AliAODTrack::kOrphan);
630 }
631 } // end of loop on tracks
632
633
634 // Access to the AOD container of vertices
635 TClonesArray &clusters = *(aod->GetNeutrals());
636 Int_t jClusters=0;
637
638 // Clusters
639 Int_t nClusters = esd->GetNumberOfCaloClusters();
640
641 for (Int_t iClust=0; iClust<nClusters; ++iClust) {
642
643 AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
644
645 Int_t id = cluster->GetID();
646 Int_t label = -1;
647 Float_t energy = cluster->GetClusterEnergy();
648 Float_t x[3];
649 cluster->GetGlobalPosition(x);
650 Float_t * covMatrix = NULL;
651 Float_t * pid = NULL;
652 AliAODVertex *prodVertex = primary;
653 AliAODTrack *primTrack = NULL;
654 Char_t ttype=AliAODNeutral::kUndef;
655
656 if (cluster->IsPHOS()) ttype=AliAODNeutral::kPHOSCluster;
657 else if (cluster->IsEMCAL()) {
658
659 if (cluster->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
660 ttype = AliAODNeutral::kEMCALPseudoCluster;
661 else
662 ttype = AliAODNeutral::kEMCALClusterv1;
663
664 }
df9db588 665
cf22b3fc 666 new(clusters[jClusters++]) AliAODNeutral(id,
667 label,
668 energy,
669 x,
670 covMatrix,
671 pid,
672 prodVertex,
673 primTrack,
674 ttype);
675
676 } // end of loop on clusters
677
678
679 delete [] usedTrack;
680 delete [] usedV0;
681 delete [] usedKink;
682
df9db588 683 // fill the tree for this event
684 aodTree->Fill();
cf22b3fc 685 } // end of event loop
df9db588 686
687 aodTree->GetUserInfo()->Add(aod);
688
689 // close ESD file
690 inFile->Close();
691
692 // write the tree to the specified file
693 outFile = aodTree->GetCurrentFile();
694 outFile->cd();
695 aodTree->Write();
696 outFile->Close();
697
698}