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