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