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