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