ANALYSIS/AliAnalysisTaskESDfilter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliAnalysisTaskESDfilter.cxx 24535 2008-03-16 22:43:30Z fca $ */
  17
  18 #include <TChain.h>
  19 #include <TFile.h>
  20 #include <TArrayI.h>
  21 #include <TRandom.h>
  22
  23 #include "AliAnalysisTaskESDfilter.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliESDEvent.h"
  26 #include "AliAODEvent.h"
  27 #include "AliESDInputHandler.h"
  28 #include "AliAODHandler.h"
  29 #include "AliAnalysisFilter.h"
  30 #include "AliESDMuonTrack.h"
  31 #include "AliESDVertex.h"
  32 #include "AliESDv0.h"
  33 #include "AliESDkink.h"
  34 #include "AliESDcascade.h"
  35 #include "AliESDPmdTrack.h"
  36 #include "AliESDCaloCluster.h"
  37 #include "AliESDCaloCells.h"
  38 #include "AliMultiplicity.h"
  39 #include "AliLog.h"
  40
  41 ClassImp(AliAnalysisTaskESDfilter)
  42
  43 ////////////////////////////////////////////////////////////////////////
  44
  45 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():
  46     AliAnalysisTaskSE(),
  47     fTrackFilter(0x0),
  48     fKinkFilter(0x0),
  49     fV0Filter(0x0),
  50     fHighPthreshold(0),
  51     fPtshape(0x0)
  52 {
  53   // Default constructor
  54 }
  55
  56 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):
  57     AliAnalysisTaskSE(name),
  58     fTrackFilter(0x0),
  59     fKinkFilter(0x0),
  60     fV0Filter(0x0),
  61     fHighPthreshold(0),
  62     fPtshape(0x0)
  63 {
  64   // Constructor
  65 }
  66
  67 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()
  68 {
  69 // Create the output container
  70     OutputTree()->GetUserInfo()->Add(fTrackFilter);
  71 }
  72
  73 void AliAnalysisTaskESDfilter::Init()
  74 {
  75     // Initialization
  76     if (fDebug > 1) AliInfo("Init() \n");
  77     // Call configuration file
  78 }
  79
  80
  81 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)
  82 {
  83 // Execute analysis for current event
  84 //
  85
  86   Long64_t ientry = Entry();
  87   if (fDebug > 0) printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
  88   if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
  89   if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
  90   ConvertESDtoAOD();
  91 }
  92
  93 void AliAnalysisTaskESDfilter::ConvertESDtoAOD() {
  94     // ESD Filter analysis task executed for each event
  95     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  96     AliESD* old = esd->GetAliESDOld();
  97
  98     // set arrays and pointers
  99     Float_t posF[3];
 100     Double_t pos[3];
 101     Double_t p[3];
 102     Double_t p_pos[3];
 103     Double_t p_neg[3];
 104     Double_t p_pos_atv0[3];
 105     Double_t p_neg_atv0[3];
 106     Double_t covVtx[6];
 107     Double_t covTr[21];
 108     Double_t pid[10];
 109
 110     for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 111     for (Int_t i = 0; i < 21; i++) covTr [i] = 0.;
 112
 113
 114   // loop over events and fill them
 115
 116   // Multiplicity information needed by the header (to be revised!)
 117     Int_t nTracks    = esd->GetNumberOfTracks();
 118     //    if (fDebug > 0) printf("-------------------Bo: Number of ESD tracks %d \n",nTracks);
 119
 120     Int_t nPosTracks = 0;
 121 //    for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack)
 122 //      if (esd->GetTrack(iTrack)->GetSign()> 0) nPosTracks++;
 123
 124     // Update the header
 125
 126     AliAODHeader* header = AODEvent()->GetHeader();
 127     header->SetRunNumber(esd->GetRunNumber());
 128     if (old) {
 129         header->SetBunchCrossNumber(0);
 130         header->SetOrbitNumber(0);
 131         header->SetPeriodNumber(0);
 132         header->SetEventType(0);
 133         header->SetMuonMagFieldScale(-999.); // FIXME
 134         header->SetCentrality(-999.);        // FIXME
 135     } else {
 136         header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
 137         header->SetOrbitNumber(esd->GetOrbitNumber());
 138         header->SetPeriodNumber(esd->GetPeriodNumber());
 139         header->SetEventType(esd->GetEventType());
 140         header->SetMuonMagFieldScale(-999.); // FIXME
 141         header->SetCentrality(-999.);        // FIXME
 142     }
 143
 144     header->SetTriggerMask(esd->GetTriggerMask());
 145     header->SetTriggerCluster(esd->GetTriggerCluster());
 146     header->SetMagneticField(esd->GetMagneticField());
 147     header->SetZDCN1Energy(esd->GetZDCN1Energy());
 148     header->SetZDCP1Energy(esd->GetZDCP1Energy());
 149     header->SetZDCN2Energy(esd->GetZDCN2Energy());
 150     header->SetZDCP2Energy(esd->GetZDCP2Energy());
 151     header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 152 //
 153 //
 154     Int_t nV0s      = esd->GetNumberOfV0s();
 155     Int_t nCascades = esd->GetNumberOfCascades();
 156     Int_t nKinks    = esd->GetNumberOfKinks();
 157     Int_t nVertices = nV0s + 2*nCascades /*could lead to two vertices, one V0 and the Xi */+ nKinks + 1 /* = prim. vtx*/;
 158     Int_t nJets     = 0;
 159     Int_t nCaloClus = esd->GetNumberOfCaloClusters();
 160     Int_t nFmdClus  = 0;
 161     Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
 162
 163     if (fDebug > 0)
 164         printf("   NV0=%d  NCASCADES=%d  NKINKS=%d\n", nV0s, nCascades, nKinks);
 165
 166     AODEvent()->ResetStd(nTracks, nVertices, nV0s+nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);
 167
 168     AliAODTrack *aodTrack = 0x0;
 169     AliAODPid   *detpid   = 0x0;
 170     Double_t timezero = 0; //TO BE FIXED
 171     AliAODv0    *aodV0    = 0x0;
 172
 173     // RefArray to take into account the tracks associated to V0s
 174     TRefArray   *v0DaughterTracks = NULL;
 175     if (nTracks>0) v0DaughterTracks = new TRefArray(nTracks);
 176
 177     // Array to take into account the tracks already added to the AOD
 178     Bool_t * usedTrack = NULL;
 179     if (nTracks>0) {
 180         usedTrack = new Bool_t[nTracks];
 181         for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) usedTrack[iTrack]=kFALSE;
 182     }
 183     // Array to take into account the V0s already added to the AOD
 184     Bool_t * usedV0 = NULL;
 185     if (nV0s>0) {
 186         usedV0 = new Bool_t[nV0s];
 187         for (Int_t iV0=0; iV0<nV0s; ++iV0) usedV0[iV0]=kFALSE;
 188     }
 189     // Array to take into account the kinks already added to the AOD
 190     Bool_t * usedKink = NULL;
 191     if (nKinks>0) {
 192         usedKink = new Bool_t[nKinks];
 193         for (Int_t iKink=0; iKink<nKinks; ++iKink) usedKink[iKink]=kFALSE;
 194     }
 195
 196     // Access to the AOD container of vertices
 197     TClonesArray &vertices = *(AODEvent()->GetVertices());
 198     Int_t jVertices=0;
 199
 200     // Access to the AOD container of tracks
 201     TClonesArray &tracks = *(AODEvent()->GetTracks());
 202     Int_t jTracks=0;
 203
 204     // Access to the AOD container of V0s
 205     TClonesArray &V0s = *(AODEvent()->GetV0s());
 206     Int_t jV0s=0;
 207
 208     // Add primary vertex. The primary tracks will be defined
 209     // after the loops on the composite objects (V0, cascades, kinks)
 210     const AliESDVertex *vtx = esd->GetPrimaryVertex();
 211
 212     vtx->GetXYZ(pos); // position
 213     vtx->GetCovMatrix(covVtx); //covariance matrix
 214
 215     AliAODVertex * primary = new(vertices[jVertices++])
 216         AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
 217     if (fDebug > 0) primary->Print();
 218
 219     // Create vertices starting from the most complex objects
 220     Double_t chi2 = 0.;
 221
 222     // Cascades
 223     for (Int_t nCascade = 0; nCascade < nCascades; ++nCascade) {
 224         AliESDcascade *cascade = esd->GetCascade(nCascade);
 225
 226         cascade->GetXYZ(pos[0], pos[1], pos[2]);
 227
 228         if (!old) {
 229             chi2 = cascade->GetChi2Xi(); // = chi2/NDF since NDF = 2*2-3
 230             cascade->GetPosCovXi(covVtx);
 231         } else {
 232             chi2 = -999.;
 233             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 234         }
 235         // Add the cascade vertex
 236         AliAODVertex * vcascade = new(vertices[jVertices++]) AliAODVertex(pos,
 237                                                                           covVtx,
 238                                                                           chi2,
 239                                                                           primary,
 240                                                                           nCascade,
 241                                                                           AliAODVertex::kCascade);
 242
 243         primary->AddDaughter(vcascade);
 244
 245         // Add the V0 from the cascade. The ESD class have to be optimized...
 246         // Now we have to search for the corresponding Vo in the list of V0s
 247         // using the indeces of the positive and negative tracks
 248
 249         Int_t posFromV0 = cascade->GetPindex();
 250         Int_t negFromV0 = cascade->GetNindex();
 251
 252
 253         AliESDv0 * v0 = 0x0;
 254         Int_t indV0 = -1;
 255
 256         for (Int_t iV0=0; iV0<nV0s; ++iV0) {
 257
 258             v0 = esd->GetV0(iV0);
 259             Int_t posV0 = v0->GetPindex();
 260             Int_t negV0 = v0->GetNindex();
 261
 262             if (posV0==posFromV0 && negV0==negFromV0) {
 263                 indV0 = iV0;
 264                 break;
 265             }
 266         }
 267
 268         AliAODVertex * vV0FromCascade = 0x0;
 269
 270         if (indV0>-1 && !usedV0[indV0] ) {
 271
 272             // the V0 exists in the array of V0s and is not used
 273
 274             usedV0[indV0] = kTRUE;
 275
 276             v0->GetXYZ(pos[0], pos[1], pos[2]);
 277             if (!old) {
 278                 chi2 = v0->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3
 279                 v0->GetPosCov(covVtx);
 280             } else {
 281                 chi2 = -999.;
 282                 for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 283             }
 284
 285             vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
 286                                                                      covVtx,
 287                                                                      chi2,
 288                                                                      vcascade,
 289                                                                      indV0,
 290                                                                      AliAODVertex::kV0);
 291         } else {
 292
 293             // the V0 doesn't exist in the array of V0s or was used
 294 //          cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
 295 //               << " The V0 " << indV0
 296 //               << " doesn't exist in the array of V0s or was used!" << endl;
 297
 298             cascade->GetXYZ(pos[0], pos[1], pos[2]);
 299
 300             if (!old) {
 301                 chi2 = v0->GetChi2V0();
 302                 cascade->GetPosCov(covVtx);
 303             } else {
 304                 chi2 = -999.;
 305                 for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 306             }
 307
 308             vV0FromCascade = new(vertices[jVertices++]) AliAODVertex(pos,
 309                                                                      covVtx,
 310                                                                      chi2, // = chi2/NDF since NDF = 2*2-3 (AM)
 311                                                                      vcascade,
 312                                                                      indV0,
 313                                                                      AliAODVertex::kV0);
 314             vcascade->AddDaughter(vV0FromCascade);
 315         }
 316
 317         // Add the positive tracks from the V0
 318
 319         if (posFromV0>-1 && !usedTrack[posFromV0]) {
 320
 321             usedTrack[posFromV0] = kTRUE;
 322
 323             AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
 324             esdTrack->GetPxPyPz(p_pos);
 325             esdTrack->GetXYZ(pos);
 326             esdTrack->GetCovarianceXYZPxPyPz(covTr);
 327             esdTrack->GetESDpid(pid);
 328
 329             vV0FromCascade->AddDaughter(aodTrack =
 330                                         new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 331                                                                            esdTrack->GetLabel(),
 332                                                                            p_pos,
 333                                                                            kTRUE,
 334                                                                            pos,
 335                                                                            kFALSE,
 336                                                                            covTr,
 337                                                                            (Short_t)esdTrack->GetSign(),
 338                                                                            esdTrack->GetITSClusterMap(),
 339                                                                            pid,
 340                                                                            vV0FromCascade,
 341                                                                            kTRUE,  // check if this is right
 342                                                                            kFALSE, // check if this is right
 343                                                                            AliAODTrack::kSecondary)
 344                 );
 345        if (esdTrack->GetSign() > 0) nPosTracks++;
 346             aodTrack->ConvertAliPIDtoAODPID();
 347             aodTrack->SetFlags(esdTrack->GetStatus());
 348             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 349         }
 350         else {
 351 //          cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
 352 //               << " track " << posFromV0 << " has already been used!" << endl;
 353         }
 354
 355         // Add the negative tracks from the V0
 356
 357         if (negFromV0>-1 && !usedTrack[negFromV0]) {
 358
 359             usedTrack[negFromV0] = kTRUE;
 360
 361             AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
 362             esdTrack->GetPxPyPz(p_neg);
 363             esdTrack->GetXYZ(pos);
 364             esdTrack->GetCovarianceXYZPxPyPz(covTr);
 365             esdTrack->GetESDpid(pid);
 366
 367             vV0FromCascade->AddDaughter(aodTrack =
 368                                         new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 369                                                                            esdTrack->GetLabel(),
 370                                                                            p_neg,
 371                                                                            kTRUE,
 372                                                                            pos,
 373                                                                            kFALSE,
 374                                                                            covTr,
 375                                                                            (Short_t)esdTrack->GetSign(),
 376                                                                            esdTrack->GetITSClusterMap(),
 377                                                                            pid,
 378                                                                            vV0FromCascade,
 379                                                                            kTRUE,  // check if this is right
 380                                                                            kFALSE, // check if this is right
 381                                                                            AliAODTrack::kSecondary)
 382                 );
 383        if (esdTrack->GetSign() > 0) nPosTracks++;
 384             aodTrack->ConvertAliPIDtoAODPID();
 385             aodTrack->SetFlags(esdTrack->GetStatus());
 386             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 387         }
 388         else {
 389 //          cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
 390 //               << " track " << negFromV0 << " has already been used!" << endl;
 391         }
 392
 393         // add it to the V0 array as well
 394         Double_t d0[2] = { -999., -99.};
 395         // counting is probably wrong
 396         new(V0s[jV0s++]) AliAODv0(vV0FromCascade, -999., -99., p_pos, p_neg, d0); // to be refined
 397
 398         // Add the bachelor track from the cascade
 399
 400         Int_t bachelor = cascade->GetBindex();
 401
 402         if(bachelor>-1 && !usedTrack[bachelor]) {
 403
 404             usedTrack[bachelor] = kTRUE;
 405
 406             AliESDtrack *esdTrack = esd->GetTrack(bachelor);
 407             esdTrack->GetPxPyPz(p);
 408             esdTrack->GetXYZ(pos);
 409             esdTrack->GetCovarianceXYZPxPyPz(covTr);
 410             esdTrack->GetESDpid(pid);
 411
 412             vcascade->AddDaughter(aodTrack =
 413                                   new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 414                                                                      esdTrack->GetLabel(),
 415                                                                      p,
 416                                                                      kTRUE,
 417                                                                      pos,
 418                                                                      kFALSE,
 419                                                                      covTr,
 420                                                                      (Short_t)esdTrack->GetSign(),
 421                                                                      esdTrack->GetITSClusterMap(),
 422                                                                      pid,
 423                                                                      vcascade,
 424                                                                      kTRUE,  // check if this is right
 425                                                                      kFALSE, // check if this is right
 426                                                                      AliAODTrack::kSecondary)
 427                 );
 428        if (esdTrack->GetSign() > 0) nPosTracks++;
 429             aodTrack->ConvertAliPIDtoAODPID();
 430             aodTrack->SetFlags(esdTrack->GetStatus());
 431             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 432         }
 433         else {
 434 //          cerr << "Error: event " << esd->GetEventNumberInFile() << " cascade " << nCascade
 435 //               << " track " << bachelor << " has already been used!" << endl;
 436         }
 437
 438         // Add the primary track of the cascade (if any)
 439
 440     } // end of the loop on cascades
 441
 442     // V0s
 443
 444     for (Int_t nV0 = 0; nV0 < nV0s; ++nV0) {
 445
 446         if (usedV0[nV0]) continue; // skip if already added to the AOD
 447
 448         AliESDv0 *v0 = esd->GetV0(nV0);
 449
 450         v0->GetXYZ(pos[0], pos[1], pos[2]);
 451
 452         if (!old) {
 453             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
 454             v0->GetPosCov(covVtx);
 455         } else {
 456             chi2 = -999.;
 457             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 458         }
 459
 460
 461         AliAODVertex * vV0 =
 462           new(vertices[jVertices++]) AliAODVertex(pos,
 463                                                   covVtx,
 464                                                   chi2,
 465                                                   primary,
 466                                                   nV0,
 467                                                   AliAODVertex::kV0);
 468         primary->AddDaughter(vV0);
 469
 470         Int_t posFromV0 = v0->GetPindex();
 471         Int_t negFromV0 = v0->GetNindex();
 472
 473         Float_t  dcaPosToPrimVertexXYZ[2]   = { 999., 999.}; // ..[0] = in XY plane and ..[1] = in Z
 474         Float_t  dcaNegToPrimVertexXYZ[2]   = { 999., 999.}; // ..[0] = in XY plane and ..[1] = in Z
 475         Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = Pos and ..[1] = Neg
 476
 477         Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();
 478         Double_t  dcaV0ToPrimVertex   = v0->GetD();
 479
 480         v0->GetPPxPyPz(p_pos_atv0[0],p_pos_atv0[1],p_pos_atv0[2]);
 481         v0->GetNPxPyPz(p_neg_atv0[0],p_neg_atv0[1],p_neg_atv0[2]);
 482
 483         // Add the positive tracks from the V0
 484
 485         if (posFromV0>-1){
 486           AliESDtrack *esdTrack = esd->GetTrack(posFromV0);
 487           esdTrack->GetPxPyPz(p_pos);
 488           esdTrack->GetXYZ(pos);
 489           esdTrack->GetCovarianceXYZPxPyPz(covTr);
 490           esdTrack->GetESDpid(pid);
 491           esdTrack->GetImpactParameters(dcaPosToPrimVertexXYZ[0],dcaPosToPrimVertexXYZ[1]);
 492           if (!usedTrack[posFromV0]) {
 493             usedTrack[posFromV0] = kTRUE;
 494             aodTrack = new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 495                                                           esdTrack->GetLabel(),
 496                                                           p_pos,
 497                                                           kTRUE,
 498                                                           pos,
 499                                                           kFALSE,
 500                                                           covTr,
 501                                                           (Short_t)esdTrack->GetSign(),
 502                                                           esdTrack->GetITSClusterMap(),
 503                                                           pid,
 504                                                           vV0,
 505                                                           kTRUE,  // check if this is right
 506                                                           kFALSE, // check if this is right
 507                                                           AliAODTrack::kSecondary);
 508             v0DaughterTracks->AddAt(aodTrack,posFromV0);
 509             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
 510             if (esdTrack->GetSign() > 0) nPosTracks++;
 511             aodTrack->ConvertAliPIDtoAODPID();
 512             aodTrack->SetFlags(esdTrack->GetStatus());
 513             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 514           }
 515           else {
 516             aodTrack = dynamic_cast<AliAODTrack*>(v0DaughterTracks->At(posFromV0));
 517             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
 518           }
 519           vV0->AddDaughter(aodTrack);
 520         }
 521
 522         // Add the negative tracks from the V0
 523
 524         if (negFromV0>-1){
 525           AliESDtrack *esdTrack = esd->GetTrack(negFromV0);
 526           esdTrack->GetPxPyPz(p_neg);
 527           esdTrack->GetXYZ(pos);
 528           esdTrack->GetCovarianceXYZPxPyPz(covTr);
 529           esdTrack->GetESDpid(pid);
 530           esdTrack->GetImpactParameters(dcaNegToPrimVertexXYZ[0],dcaNegToPrimVertexXYZ[1]);
 531
 532           if (!usedTrack[negFromV0]) {
 533             usedTrack[negFromV0] = kTRUE;
 534             aodTrack = new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 535                                                           esdTrack->GetLabel(),
 536                                                           p_neg,
 537                                                           kTRUE,
 538                                                           pos,
 539                                                           kFALSE,
 540                                                           covTr,
 541                                                           (Short_t)esdTrack->GetSign(),
 542                                                           esdTrack->GetITSClusterMap(),
 543                                                           pid,
 544                                                           vV0,
 545                                                           kTRUE,  // check if this is right
 546                                                           kFALSE, // check if this is right
 547                                                           AliAODTrack::kSecondary);
 548
 549             v0DaughterTracks->AddAt(aodTrack,negFromV0);
 550             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
 551             if (esdTrack->GetSign() > 0) nPosTracks++;
 552             aodTrack->ConvertAliPIDtoAODPID();
 553             aodTrack->SetFlags(esdTrack->GetStatus());
 554             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 555           }
 556           else {
 557             aodTrack = dynamic_cast<AliAODTrack*>(v0DaughterTracks->At(negFromV0));
 558             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
 559           }
 560           vV0->AddDaughter(aodTrack);
 561         }
 562         dcaDaughterToPrimVertex[0] =
 563           TMath::Sqrt(dcaPosToPrimVertexXYZ[0]*dcaPosToPrimVertexXYZ[0]
 564                       +dcaPosToPrimVertexXYZ[1]*dcaPosToPrimVertexXYZ[1]);
 565         dcaDaughterToPrimVertex[1] =
 566           TMath::Sqrt(dcaNegToPrimVertexXYZ[0]*dcaNegToPrimVertexXYZ[0]
 567                       +dcaNegToPrimVertexXYZ[1]*dcaNegToPrimVertexXYZ[1]);
 568         // add it to the V0 array as well
 569         aodV0 = new(V0s[jV0s++])
 570           AliAODv0(vV0, dcaV0Daughters, dcaV0ToPrimVertex, p_pos_atv0, p_neg_atv0, dcaDaughterToPrimVertex); // to be refined
 571         // set the aod v0 on-the-fly status
 572         aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
 573 //      if (fDebug > 0)
 574 //        printf("-------------------Bo: In the v0 loop: onFlyStatus=%d,dcaV0Daughters=%.3f dcaV0ToPrimVertex=%.3f dcaDaughterToPrimVertex=[%.3f,%.3f]\n",
 575 //               aodV0->GetOnFlyStatus(), aodV0->DcaV0Daughters(), aodV0->DcaV0ToPrimVertex(),
 576 //               aodV0->DcaPosToPrimVertex(),aodV0->DcaNegToPrimVertex());
 577     }
 578     V0s.Expand(jV0s);
 579     // end of the loop on V0s
 580
 581     // Kinks: it is a big mess the access to the information in the kinks
 582     // The loop is on the tracks in order to find the mother and daugther of each kink
 583
 584
 585     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
 586
 587         AliESDtrack * esdTrack = esd->GetTrack(iTrack);
 588
 589         Int_t ikink = esdTrack->GetKinkIndex(0);
 590
 591         if (ikink && nKinks) {
 592             // Negative kink index: mother, positive: daughter
 593
 594             // Search for the second track of the kink
 595
 596             for (Int_t jTrack = iTrack+1; jTrack<nTracks; ++jTrack) {
 597
 598                 AliESDtrack * esdTrack1 = esd->GetTrack(jTrack);
 599
 600                 Int_t jkink = esdTrack1->GetKinkIndex(0);
 601
 602                 if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
 603
 604                     // The two tracks are from the same kink
 605
 606                     if (usedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
 607
 608                     Int_t imother = -1;
 609                     Int_t idaughter = -1;
 610
 611                     if (ikink<0 && jkink>0) {
 612
 613                         imother = iTrack;
 614                         idaughter = jTrack;
 615                     }
 616                     else if (ikink>0 && jkink<0) {
 617
 618                         imother = jTrack;
 619                         idaughter = iTrack;
 620                     }
 621                     else {
 622 //                      cerr << "Error: Wrong combination of kink indexes: "
 623 //                           << ikink << " " << jkink << endl;
 624                         continue;
 625                     }
 626
 627                     // Add the mother track
 628
 629                     AliAODTrack * mother = NULL;
 630
 631                     if (!usedTrack[imother]) {
 632
 633                         usedTrack[imother] = kTRUE;
 634
 635                         AliESDtrack *esdTrackM = esd->GetTrack(imother);
 636                         esdTrackM->GetPxPyPz(p);
 637                         esdTrackM->GetXYZ(pos);
 638                         esdTrackM->GetCovarianceXYZPxPyPz(covTr);
 639                         esdTrackM->GetESDpid(pid);
 640
 641                         mother =
 642                             new(tracks[jTracks++]) AliAODTrack(esdTrackM->GetID(),
 643                                                                esdTrackM->GetLabel(),
 644                                                                p,
 645                                                                kTRUE,
 646                                                                pos,
 647                                                                kFALSE,
 648                                                                covTr,
 649                                                                (Short_t)esdTrackM->GetSign(),
 650                                                                esdTrackM->GetITSClusterMap(),
 651                                                                pid,
 652                                                                primary,
 653                                                                kTRUE, // check if this is right
 654                                                                kTRUE, // check if this is right
 655                                                                AliAODTrack::kPrimary);
 656                         if (esdTrackM->GetSign() > 0) nPosTracks++;
 657                         mother->SetFlags(esdTrackM->GetStatus());
 658                         mother->ConvertAliPIDtoAODPID();
 659                         primary->AddDaughter(mother);
 660                         mother->ConvertAliPIDtoAODPID();
 661                         SetAODPID(esdTrackM,mother,detpid,timezero);
 662                     }
 663                     else {
 664 //                      cerr << "Error: event " << esd->GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
 665 //                           << " track " << imother << " has already been used!" << endl;
 666                   }
 667
 668                     // Add the kink vertex
 669                     AliESDkink * kink = esd->GetKink(TMath::Abs(ikink)-1);
 670
 671                     AliAODVertex * vkink =
 672                         new(vertices[jVertices++]) AliAODVertex(kink->GetPosition(),
 673                                                                 NULL,
 674                                                                 0.,
 675                                                                 mother,
 676                                                                 esdTrack->GetID(),  // This is the track ID of the mother's track!
 677                                                                 AliAODVertex::kKink);
 678                     // Add the daughter track
 679
 680                     AliAODTrack * daughter = NULL;
 681
 682                     if (!usedTrack[idaughter]) {
 683
 684                         usedTrack[idaughter] = kTRUE;
 685
 686                         AliESDtrack *esdTrackD = esd->GetTrack(idaughter);
 687                         esdTrackD->GetPxPyPz(p);
 688                         esdTrackD->GetXYZ(pos);
 689                         esdTrackD->GetCovarianceXYZPxPyPz(covTr);
 690                         esdTrackD->GetESDpid(pid);
 691
 692                         daughter =
 693                             new(tracks[jTracks++]) AliAODTrack(esdTrackD->GetID(),
 694                                                                esdTrackD->GetLabel(),
 695                                                                p,
 696                                                                kTRUE,
 697                                                                pos,
 698                                                                kFALSE,
 699                                                                covTr,
 700                                                                (Short_t)esdTrackD->GetSign(),
 701                                                                esdTrackD->GetITSClusterMap(),
 702                                                                pid,
 703                                                                vkink,
 704                                                                kTRUE, // check if this is right
 705                                                                kTRUE, // check if this is right
 706                                                                AliAODTrack::kSecondary);
 707                         if (esdTrackD->GetSign() > 0) nPosTracks++;
 708                         daughter->SetFlags(esdTrackD->GetStatus());
 709                         daughter->ConvertAliPIDtoAODPID();
 710                         vkink->AddDaughter(daughter);
 711                         daughter->ConvertAliPIDtoAODPID();
 712                        SetAODPID(esdTrackD,daughter,detpid,timezero);
 713                     }
 714                     else {
 715 //                      cerr << "Error: event " << esd->GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
 716 //                           << " track " << idaughter << " has already been used!" << endl;
 717                     }
 718                 }
 719             }
 720         }
 721     }
 722
 723
 724     // Tracks (primary and orphan)
 725
 726     if (fDebug > 0) printf("NUMBER OF ESD TRACKS %5d\n", nTracks);
 727
 728     for (Int_t nTrack = 0; nTrack < nTracks; ++nTrack) {
 729
 730
 731         if (usedTrack[nTrack]) continue;
 732
 733         AliESDtrack *esdTrack = esd->GetTrack(nTrack);
 734         UInt_t selectInfo = 0;
 735         //
 736         // Track selection
 737         if (fTrackFilter) {
 738             selectInfo = fTrackFilter->IsSelected(esdTrack);
 739             if (!selectInfo) continue;
 740         }
 741
 742         //
 743         esdTrack->GetPxPyPz(p);
 744         esdTrack->GetXYZ(pos);
 745         esdTrack->GetCovarianceXYZPxPyPz(covTr);
 746         esdTrack->GetESDpid(pid);
 747
 748         Float_t impactXY, impactZ;
 749
 750         esdTrack->GetImpactParameters(impactXY,impactZ);
 751
 752         if (impactXY<3) {
 753             // track inside the beam pipe
 754
 755             primary->AddDaughter(aodTrack =
 756                                  new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 757                                                                     esdTrack->GetLabel(),
 758                                                                     p,
 759                                                                     kTRUE,
 760                                                                     pos,
 761                                                                     kFALSE,
 762                                                                     covTr,
 763                                                                     (Short_t)esdTrack->GetSign(),
 764                                                                     esdTrack->GetITSClusterMap(),
 765                                                                     pid,
 766                                                                     primary,
 767                                                                     kTRUE, // check if this is right
 768                                                                     kTRUE, // check if this is right
 769                                                                     AliAODTrack::kPrimary,
 770                                                                     selectInfo)
 771                 );
 772        if (esdTrack->GetSign() > 0) nPosTracks++;
 773             aodTrack->SetFlags(esdTrack->GetStatus());
 774             aodTrack->ConvertAliPIDtoAODPID();
 775             SetAODPID(esdTrack,aodTrack,detpid,timezero);
 776         }
 777         else {
 778           // outside the beam pipe: orphan track
 779           // Don't write them anymore!
 780           continue;
 781           /*
 782           aodTrack =
 783             new(tracks[jTracks++]) AliAODTrack(esdTrack->GetID(),
 784                                                esdTrack->GetLabel(),
 785                                                p,
 786                                                kTRUE,
 787                                                pos,
 788                                                kFALSE,
 789                                                covTr,
 790                                                (Short_t)esdTrack->GetSign(),
 791                                                esdTrack->GetITSClusterMap(),
 792                                                pid,
 793                                                NULL,
 794                                                kFALSE, // check if this is right
 795                                                kFALSE, // check if this is right
 796                                                AliAODTrack::kOrphan,
 797                                                selectInfo);
 798      if (esdTrack->GetSign() > 0) nPosTracks++;
 799           aodTrack->SetFlags(esdTrack->GetStatus());
 800           aodTrack->ConvertAliPIDtoAODPID();
 801           */
 802         }
 803     } // end of loop on tracks
 804
 805     // Update number of AOD tracks in header at the end of track loop (M.G.)
 806     header->SetRefMultiplicity(jTracks);
 807     header->SetRefMultiplicityPos(nPosTracks);
 808     header->SetRefMultiplicityNeg(jTracks - nPosTracks);
 809     if (fDebug > 0)
 810       printf("   NAODTRACKS=%d  NPOS=%d  NNEG=%d\n", jTracks, nPosTracks, jTracks - nPosTracks);
 811     // Do not shrink the array of tracks - other filters may add to it (M.G)
 812 //    tracks.Expand(jTracks); // remove 'empty slots' due to unwritten tracks
 813
 814     // Access to the AOD container of PMD clusters
 815     TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
 816     Int_t jPmdClusters=0;
 817
 818     for (Int_t iPmd = 0; iPmd < nPmdClus; ++iPmd) {
 819       // file pmd clusters, to be revised!
 820       AliESDPmdTrack *pmdTrack = esd->GetPmdTrack(iPmd);
 821       Int_t nLabel = 0;
 822       Int_t *label = 0x0;
 823       Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
 824       Double_t pidPmd[9] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; // to be revised!
 825       // type not set!
 826       // assoc cluster not set
 827       new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
 828     }
 829
 830     // Access to the AOD container of clusters
 831     TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 832     Int_t jClusters=0;
 833
 834     for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 835
 836       AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
 837
 838       Int_t id        = cluster->GetID();
 839       Int_t nLabel    = cluster->GetNLabels();
 840       TArrayI* labels = cluster->GetLabels();
 841       Int_t *label = 0;
 842       if (labels) label = (cluster->GetLabels())->GetArray();
 843
 844       Float_t energy = cluster->E();
 845       cluster->GetPosition(posF);
 846       Char_t ttype = AliAODCluster::kUndef;
 847
 848       if (cluster->GetClusterType() == AliESDCaloCluster::kPHOSCluster) {
 849         ttype=AliAODCluster::kPHOSNeutral;
 850       }
 851       else if (cluster->GetClusterType() == AliESDCaloCluster::kEMCALClusterv1) {
 852         ttype = AliAODCluster::kEMCALClusterv1;
 853       }
 854
 855
 856       AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
 857                                                                                         nLabel,
 858                                                                                         label,
 859                                                                                         energy,
 860                                                                                         pos,
 861                                                                                         NULL,
 862                                                                                         ttype);
 863
 864       caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 865                                   cluster->GetClusterDisp(),
 866                                   cluster->GetM20(), cluster->GetM02(),
 867                                   cluster->GetEmcCpvDistance(),
 868                                   cluster->GetNExMax(),cluster->GetTOF()) ;
 869
 870       caloCluster->SetPIDFromESD(cluster->GetPid());
 871       caloCluster->SetNCells(cluster->GetNCells());
 872       caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 873       caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 874
 875       TArrayI* matchedT =       cluster->GetTracksMatched();
 876       if (matchedT) {
 877         for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 878           caloCluster->AddTrackMatched((esd->GetTrack(im)));
 879         }
 880       }
 881
 882     }
 883     caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 884     // end of loop on calo clusters
 885
 886     // fill EMCAL cell info
 887     if (esd->GetEMCALCells()) { // protection against missing ESD information
 888       AliESDCaloCells &esdEMcells = *(esd->GetEMCALCells());
 889       Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
 890
 891       AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 892       aodEMcells.CreateContainer(nEMcell);
 893       aodEMcells.SetType(AliAODCaloCells::kEMCAL);
 894       for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 895         aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell));
 896       }
 897       aodEMcells.Sort();
 898     }
 899
 900     // fill PHOS cell info
 901     if (esd->GetPHOSCells()) { // protection against missing ESD information
 902       AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
 903       Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
 904
 905       AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 906       aodPHcells.CreateContainer(nPHcell);
 907       aodPHcells.SetType(AliAODCaloCells::kPHOS);
 908       for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 909         aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
 910       }
 911       aodPHcells.Sort();
 912     }
 913
 914     // tracklets
 915     AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());
 916     const AliMultiplicity *mult = esd->GetMultiplicity();
 917     if (mult) {
 918       if (mult->GetNumberOfTracklets()>0) {
 919         SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
 920
 921         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
 922           SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0), mult->GetLabel(n, 1));
 923         }
 924       }
 925     } else {
 926       //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
 927     }
 928
 929     delete [] usedTrack;
 930     delete [] usedV0;
 931     delete [] usedKink;
 932     delete    v0DaughterTracks;
 933
 934     return;
 935 }
 936
 937 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid, Double_t timezero)
 938 {
 939   //
 940   // Setter for the raw PID detector signals
 941   //
 942
 943   if(esdtrack->Pt()>fHighPthreshold) {
 944     detpid = new AliAODPid();
 945     detpid->SetDetectorRawSignals(esdtrack,timezero);
 946     aodtrack->SetDetPID(detpid);
 947   } else {
 948     if(fPtshape){
 949       if(esdtrack->Pt()> fPtshape->GetXmin()){
 950         Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
 951         if(gRandom->Rndm(0)<1./y){
 952           detpid = new AliAODPid();
 953           detpid->SetDetectorRawSignals(esdtrack,timezero);
 954           aodtrack->SetDetPID(detpid);
 955         }//end rndm
 956       }//end if p < pmin
 957     }//end if p function
 958   }// end else
 959 }
 960
 961
 962 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
 963 {
 964 // Terminate analysis
 965 //
 966     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
 967 }
 968