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$ */
  17
  18 #include <TChain.h>
  19 #include <TTree.h>
  20 #include <TList.h>
  21 #include <TArrayI.h>
  22 #include <TParameter.h>
  23 #include <TRandom.h>
  24 #include <TParticle.h>
  25 #include <TFile.h>
  26
  27 #include "AliAnalysisTaskESDfilter.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliESDEvent.h"
  30 #include "AliESDRun.h"
  31 #include "AliStack.h"
  32 #include "AliAODEvent.h"
  33 #include "AliMCEvent.h"
  34 #include "AliMCEventHandler.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAODMCParticle.h"
  38 #include "AliAnalysisFilter.h"
  39 #include "AliESDMuonTrack.h"
  40 #include "AliESDVertex.h"
  41 #include "AliCentrality.h"
  42 #include "AliEventplane.h"
  43 #include "AliESDv0.h"
  44 #include "AliESDkink.h"
  45 #include "AliESDcascade.h"
  46 #include "AliESDPmdTrack.h"
  47 #include "AliESDCaloCluster.h"
  48 #include "AliESDCaloCells.h"
  49 #include "AliMultiplicity.h"
  50 #include "AliLog.h"
  51 #include "AliCodeTimer.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "AliESDpid.h"
  54 #include "AliAODHMPIDrings.h"
  55 #include "AliV0vertexer.h"
  56 #include "AliCascadeVertexer.h"
  57 #include "Riostream.h"
  58 #include "AliExternalTrackParam.h"
  59 #include "AliTrackerBase.h"
  60 #include "TVector3.h"
  61 #include "AliTPCdEdxInfo.h"
  62
  63 using std::cout;
  64 using std::endl;
  65 ClassImp(AliAnalysisTaskESDfilter)
  66
  67 ////////////////////////////////////////////////////////////////////////
  68
  69 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():
  70   AliAnalysisTaskSE(),
  71   fTrackFilter(0x0),
  72   fKinkFilter(0x0),
  73   fV0Filter(0x0),
  74   fCascadeFilter(0x0),
  75   fHighPthreshold(0),
  76   fPtshape(0x0),
  77   fEnableFillAOD(kTRUE),
  78   fUsedTrack(0x0),
  79   fUsedKink(0x0),
  80   fUsedV0(0x0),
  81   fAODTrackRefs(0x0),
  82   fAODV0VtxRefs(0x0),
  83   fAODV0Refs(0x0),
  84   fMChandler(0x0),
  85   fNumberOfTracks(0),
  86   fNumberOfPositiveTracks(0),
  87   fNumberOfV0s(0),
  88   fNumberOfVertices(0),
  89   fNumberOfCascades(0),
  90   fNumberOfKinks(0),
  91   fOldESDformat(kFALSE),
  92   fPrimaryVertex(0x0),
  93   fTPCConstrainedFilterMask(0),
  94   fHybridFilterMaskTPCCG(0),
  95   fWriteHybridTPCCOnly(kFALSE),
  96   fGlobalConstrainedFilterMask(0),
  97   fHybridFilterMaskGCG(0),
  98   fWriteHybridGCOnly(kFALSE),
  99   fIsVZEROEnabled(kTRUE),
 100   fIsTZEROEnabled(kTRUE),
 101   fIsZDCEnabled(kTRUE),
 102   fIsHMPIDEnabled(kTRUE),
 103   fIsV0CascadeRecoEnabled(kFALSE),
 104   fAreCascadesEnabled(kTRUE),
 105   fAreV0sEnabled(kTRUE),
 106   fAreKinksEnabled(kTRUE),
 107   fAreTracksEnabled(kTRUE),
 108   fArePmdClustersEnabled(kTRUE),
 109   fAreCaloClustersEnabled(kTRUE),
 110   fAreEMCALCellsEnabled(kTRUE),
 111   fArePHOSCellsEnabled(kTRUE),
 112   fAreEMCALTriggerEnabled(kTRUE),
 113   fArePHOSTriggerEnabled(kTRUE),
 114   fAreTrackletsEnabled(kTRUE),
 115   fESDpid(0x0),
 116   fIsPidOwner(kFALSE),
 117   fTPCaloneTrackCuts(0),
 118   fDoPropagateTrackToEMCal(kTRUE)
 119 {
 120   // Default constructor
 121     fV0Cuts[0] =  33.   ;   // max allowed chi2
 122     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 123     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 124     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 125     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 126     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 127     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 128
 129     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 130     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 131     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 132     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 133     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 134     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 135     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 136     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 137 }
 138
 139 //______________________________________________________________________________
 140 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):
 141     AliAnalysisTaskSE(name),
 142     fTrackFilter(0x0),
 143     fKinkFilter(0x0),
 144     fV0Filter(0x0),
 145     fCascadeFilter(0x0),
 146     fHighPthreshold(0),
 147     fPtshape(0x0),
 148     fEnableFillAOD(kTRUE),
 149     fUsedTrack(0x0),
 150     fUsedKink(0x0),
 151     fUsedV0(0x0),
 152     fAODTrackRefs(0x0),
 153     fAODV0VtxRefs(0x0),
 154     fAODV0Refs(0x0),
 155     fMChandler(0x0),
 156     fNumberOfTracks(0),
 157     fNumberOfPositiveTracks(0),
 158     fNumberOfV0s(0),
 159     fNumberOfVertices(0),
 160     fNumberOfCascades(0),
 161     fNumberOfKinks(0),
 162     fOldESDformat(kFALSE),
 163     fPrimaryVertex(0x0),
 164   fTPCConstrainedFilterMask(0),
 165   fHybridFilterMaskTPCCG(0),
 166   fWriteHybridTPCCOnly(kFALSE),
 167   fGlobalConstrainedFilterMask(0),
 168   fHybridFilterMaskGCG(0),
 169   fWriteHybridGCOnly(kFALSE),
 170     fIsVZEROEnabled(kTRUE),
 171     fIsTZEROEnabled(kTRUE),
 172     fIsZDCEnabled(kTRUE),
 173     fIsHMPIDEnabled(kTRUE),
 174     fIsV0CascadeRecoEnabled(kFALSE),
 175     fAreCascadesEnabled(kTRUE),
 176     fAreV0sEnabled(kTRUE),
 177     fAreKinksEnabled(kTRUE),
 178     fAreTracksEnabled(kTRUE),
 179     fArePmdClustersEnabled(kTRUE),
 180     fAreCaloClustersEnabled(kTRUE),
 181     fAreEMCALCellsEnabled(kTRUE),
 182     fArePHOSCellsEnabled(kTRUE),
 183                 fAreEMCALTriggerEnabled(kTRUE),
 184                 fArePHOSTriggerEnabled(kTRUE),
 185                 fAreTrackletsEnabled(kTRUE),
 186     fESDpid(0x0),
 187     fIsPidOwner(kFALSE),
 188     fTPCaloneTrackCuts(0),
 189   fDoPropagateTrackToEMCal(kTRUE)
 190 {
 191   // Constructor
 192
 193     fV0Cuts[0] =  33.   ;   // max allowed chi2
 194     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 195     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 196     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 197     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 198     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 199     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 200
 201     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 202     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 203     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 204     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 205     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 206     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 207     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 208     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 209
 210
 211
 212 }
 213 AliAnalysisTaskESDfilter::~AliAnalysisTaskESDfilter(){
 214     if(fIsPidOwner) delete fESDpid;
 215 }
 216 //______________________________________________________________________________
 217 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()
 218 {
 219   //
 220   // Create Output Objects conenct filter to outputtree
 221   //
 222   if(OutputTree())
 223   {
 224     OutputTree()->GetUserInfo()->Add(fTrackFilter);
 225   }
 226   else
 227   {
 228     AliError("No OutputTree() for adding the track filter");
 229   }
 230   fTPCaloneTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 231 }
 232
 233 //______________________________________________________________________________
 234 void AliAnalysisTaskESDfilter::Init()
 235 {
 236   // Initialization
 237   if (fDebug > 1) AliInfo("Init() \n");
 238   // Call configuration file
 239 }
 240
 241 //______________________________________________________________________________
 242 Bool_t AliAnalysisTaskESDfilter::Notify()
 243 {
 244 // Notify method.
 245    AddMetadataToUserInfo();
 246    return kTRUE;
 247 }
 248
 249 //______________________________________________________________________________
 250 Bool_t AliAnalysisTaskESDfilter::AddMetadataToUserInfo()
 251 {
 252 // Copy metadata to AOD user info.
 253   static Bool_t copyFirst = kFALSE;
 254   if (!copyFirst) {
 255     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 256     if (!mgr) {
 257        AliError("AliAnalysisTaskESDfilter::AddMetadataToUserInfo() : No analysis manager !");
 258        return kFALSE;
 259     }
 260     TTree *esdTree = mgr->GetTree()->GetTree();
 261     if (!esdTree) return kFALSE;
 262     TNamed *alirootVersion = (TNamed*)esdTree->GetUserInfo()->FindObject("alirootVersion");
 263     if (!alirootVersion) return kFALSE;
 264     AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
 265          if (!aodHandler) return kFALSE;
 266     TTree *aodTree = aodHandler->GetTree();
 267     if (!aodTree) return kFALSE;
 268     aodTree->GetUserInfo()->Add(new TNamed(*alirootVersion));
 269     copyFirst = kTRUE;
 270   }
 271   return kTRUE;
 272 }
 273
 274 //______________________________________________________________________________
 275 void AliAnalysisTaskESDfilter::PrintTask(Option_t *option, Int_t indent) const
 276 {
 277 // Print selection task information
 278   AliInfo("");
 279
 280   AliAnalysisTaskSE::PrintTask(option,indent);
 281
 282   TString spaces(' ',indent+3);
 283
 284         cout << spaces.Data() << Form("Cascades       are %s",fAreCascadesEnabled ? "ENABLED":"DISABLED") << endl;
 285         cout << spaces.Data() << Form("V0s            are %s",fAreV0sEnabled ? "ENABLED":"DISABLED") << endl;
 286         cout << spaces.Data() << Form("Kinks          are %s",fAreKinksEnabled ? "ENABLED":"DISABLED") << endl;
 287         cout << spaces.Data() << Form("Tracks         are %s",fAreTracksEnabled ? "ENABLED":"DISABLED") << endl;
 288         cout << spaces.Data() << Form("PmdClusters    are %s",fArePmdClustersEnabled ? "ENABLED":"DISABLED") << endl;
 289         cout << spaces.Data() << Form("CaloClusters   are %s",fAreCaloClustersEnabled ? "ENABLED":"DISABLED") << endl;
 290   cout << spaces.Data() << Form("EMCAL cells    are %s",fAreEMCALCellsEnabled ? "ENABLED":"DISABLED") << endl;
 291         cout << spaces.Data() << Form("EMCAL triggers are %s",fAreEMCALTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 292         cout << spaces.Data() << Form("PHOS triggers  are %s",fArePHOSTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 293         cout << spaces.Data() << Form("Tracklets      are %s",fAreTrackletsEnabled ? "ENABLED":"DISABLED") << endl;
 294         cout << spaces.Data() << Form("PropagateTrackToEMCal  is %s", fDoPropagateTrackToEMCal ? "ENABLED":"DISABLED") << endl;
 295 }
 296
 297 //______________________________________________________________________________
 298 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)
 299 {
 300 // Execute analysis for current event
 301 //
 302
 303   Long64_t ientry = Entry();
 304
 305   if (fDebug > 0) {
 306       printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
 307       if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
 308       if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
 309   }
 310   // Filters must explicitely enable AOD filling in their UserExec (AG)
 311   if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
 312   if(fEnableFillAOD) {
 313      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 314      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 315   }
 316   ConvertESDtoAOD();
 317 }
 318
 319 //______________________________________________________________________________
 320 TClonesArray& AliAnalysisTaskESDfilter::Cascades()
 321 {
 322   return *(AODEvent()->GetCascades());
 323 }
 324
 325 //______________________________________________________________________________
 326 TClonesArray& AliAnalysisTaskESDfilter::Tracks()
 327 {
 328   return *(AODEvent()->GetTracks());
 329 }
 330
 331 //______________________________________________________________________________
 332 TClonesArray& AliAnalysisTaskESDfilter::V0s()
 333 {
 334   return *(AODEvent()->GetV0s());
 335 }
 336
 337 //______________________________________________________________________________
 338 TClonesArray& AliAnalysisTaskESDfilter::Vertices()
 339 {
 340   return *(AODEvent()->GetVertices());
 341 }
 342
 343 //______________________________________________________________________________
 344 AliAODHeader* AliAnalysisTaskESDfilter::ConvertHeader(const AliESDEvent& esd)
 345 {
 346 // Convert header information
 347
 348   AliCodeTimerAuto("",0);
 349
 350   AliAODHeader* header = AODEvent()->GetHeader();
 351
 352   header->SetRunNumber(esd.GetRunNumber());
 353   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 354
 355   TTree* tree = fInputHandler->GetTree();
 356   if (tree) {
 357     TFile* file = tree->GetCurrentFile();
 358     if (file) header->SetESDFileName(file->GetName());
 359   }
 360
 361   if (fOldESDformat) {
 362     header->SetBunchCrossNumber(0);
 363     header->SetOrbitNumber(0);
 364     header->SetPeriodNumber(0);
 365     header->SetEventType(0);
 366     header->SetMuonMagFieldScale(-999.);
 367     header->SetCentrality(0);
 368     header->SetEventplane(0);
 369   } else {
 370     header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
 371     header->SetOrbitNumber(esd.GetOrbitNumber());
 372     header->SetPeriodNumber(esd.GetPeriodNumber());
 373     header->SetEventType(esd.GetEventType());
 374
 375     header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
 376     if(const_cast<AliESDEvent&>(esd).GetCentrality()){
 377       header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
 378     }
 379     else{
 380       header->SetCentrality(0);
 381     }
 382     if(const_cast<AliESDEvent&>(esd).GetEventplane()){
 383       header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
 384     }
 385     else{
 386       header->SetEventplane(0);
 387     }
 388   }
 389
 390   // Trigger
 391   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
 392   header->SetTriggerMask(esd.GetTriggerMask());
 393   header->SetTriggerCluster(esd.GetTriggerCluster());
 394   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());
 395   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());
 396   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());
 397
 398   header->SetMagneticField(esd.GetMagneticField());
 399   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
 400   header->SetZDCN1Energy(esd.GetZDCN1Energy());
 401   header->SetZDCP1Energy(esd.GetZDCP1Energy());
 402   header->SetZDCN2Energy(esd.GetZDCN2Energy());
 403   header->SetZDCP2Energy(esd.GetZDCP2Energy());
 404   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
 405
 406   header->SetIRInt2InteractionMap(esd.GetHeader()->GetIRInt2InteractionMap());
 407   header->SetIRInt1InteractionMap(esd.GetHeader()->GetIRInt1InteractionMap());
 408
 409   // ITS Cluster Multiplicty
 410   const AliMultiplicity *mult = esd.GetMultiplicity();
 411   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
 412
 413   // TPC only Reference Multiplicty
 414   Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
 415   header->SetTPConlyRefMultiplicity(refMult);
 416   //
 417   header->SetRefMultiplicityComb05(AliESDtrackCuts::GetReferenceMultiplicity(&esd,AliESDtrackCuts::kTrackletsITSTPC,0.5));
 418   header->SetRefMultiplicityComb08(AliESDtrackCuts::GetReferenceMultiplicity(&esd,AliESDtrackCuts::kTrackletsITSTPC,0.8));
 419   //
 420   Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};
 421   Float_t diamcov[3];
 422   esd.GetDiamondCovXY(diamcov);
 423   header->SetDiamond(diamxy,diamcov);
 424   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
 425
 426   // VZERO channel equalization factors for event-plane reconstruction
 427   header->SetVZEROEqFactors(esd.GetVZEROEqFactors());
 428
 429   // T0 Resolution information
 430   const AliESDRun* esdRun = esd.GetESDRun();
 431   for (Int_t i=0;i<AliESDRun::kT0spreadSize;i++) header->SetT0spread(i,esdRun->GetT0spread(i));
 432
 433   return header;
 434 }
 435
 436 //______________________________________________________________________________
 437 void AliAnalysisTaskESDfilter::ConvertCascades(const AliESDEvent& esd)
 438 {
 439
 440   // Convert the cascades part of the ESD.
 441   // Return the number of cascades
 442
 443   AliCodeTimerAuto("",0);
 444
 445   // Create vertices starting from the most complex objects
 446   Double_t chi2 = 0.;
 447
 448   const AliESDVertex* vtx = esd.GetPrimaryVertex();
 449   Double_t pos[3] = { 0. };
 450   Double_t covVtx[6] = { 0. };
 451   Double_t momBach[3]={0.};
 452   Double_t covTr[21]={0.};
 453   Double_t pid[10]={0.};
 454   AliAODPid* detpid(0x0);
 455   AliAODVertex* vV0FromCascade(0x0);
 456   AliAODv0* aodV0(0x0);
 457   AliAODcascade* aodCascade(0x0);
 458   AliAODTrack* aodTrack(0x0);
 459   Double_t momPos[3]={0.};
 460   Double_t momNeg[3] = { 0. };
 461   Double_t momPosAtV0vtx[3]={0.};
 462   Double_t momNegAtV0vtx[3]={0.};
 463
 464   TClonesArray& verticesArray = Vertices();
 465   TClonesArray& tracksArray = Tracks();
 466   TClonesArray& cascadesArray = Cascades();
 467
 468   // Cascades (Modified by A.Maire - February 2009)
 469   for (Int_t nCascade = 0; nCascade < esd.GetNumberOfCascades(); ++nCascade) {
 470
 471     // 0- Preparation
 472     //
 473     AliESDcascade *esdCascade = esd.GetCascade(nCascade);
 474                 Int_t  idxPosFromV0Dghter  = esdCascade->GetPindex();
 475                 Int_t  idxNegFromV0Dghter  = esdCascade->GetNindex();
 476                 Int_t  idxBachFromCascade  = esdCascade->GetBindex();
 477
 478     AliESDtrack  *esdCascadePos  = esd.GetTrack( idxPosFromV0Dghter);
 479     AliESDtrack  *esdCascadeNeg  = esd.GetTrack( idxNegFromV0Dghter);
 480     AliESDtrack  *esdCascadeBach = esd.GetTrack( idxBachFromCascade);
 481
 482     // Identification of the V0 within the esdCascade (via both daughter track indices)
 483     AliESDv0 * currentV0   = 0x0;
 484     Int_t      idxV0FromCascade = -1;
 485
 486     for (Int_t iV0=0; iV0<esd.GetNumberOfV0s(); ++iV0) {
 487
 488       currentV0 = esd.GetV0(iV0);
 489       Int_t posCurrentV0 = currentV0->GetPindex();
 490       Int_t negCurrentV0 = currentV0->GetNindex();
 491
 492       if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {
 493         idxV0FromCascade = iV0;
 494         break;
 495       }
 496     }
 497
 498     if(idxV0FromCascade < 0){
 499       printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");
 500       continue;
 501     }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"
 502
 503     AliESDv0 *esdV0FromCascade   = esd.GetV0(idxV0FromCascade);
 504
 505     // 1 - Cascade selection
 506
 507     //  AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd.GetPrimaryVertex()));
 508     //  TList cascadeObjects;
 509     //  cascadeObjects.AddAt(esdV0FromCascade, 0);
 510     //  cascadeObjects.AddAt(esdCascadePos,    1);
 511     //  cascadeObjects.AddAt(esdCascadeNeg,    2);
 512     //  cascadeObjects.AddAt(esdCascade,       3);
 513     //  cascadeObjects.AddAt(esdCascadeBach,   4);
 514     //  cascadeObjects.AddAt(esdPrimVtx,       5);
 515     //
 516     //  UInt_t selectCascade = 0;
 517     //  if (fCascadeFilter) {
 518     //    // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects);
 519     //          // FIXME AliESDCascadeCuts to be implemented ...
 520     //
 521     //          // Here we may encounter a moot point at the V0 level
 522     //          // between the cascade selections and the V0 ones :
 523     //          // the V0 selected along with the cascade (secondary V0) may
 524     //          // usually be removed from the dedicated V0 selections (prim V0) ...
 525     //          // -> To be discussed !
 526     //
 527     //    // this is a little awkward but otherwise the
 528     //    // list wants to access the pointer (delete it)
 529     //    // again when going out of scope
 530     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 531     //    esdPrimVtx = 0;
 532     //    if (!selectCascade)
 533     //      continue;
 534     //  }
 535     //  else{
 536     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 537     //    esdPrimVtx = 0;
 538     //  }
 539
 540     // 2 - Add the cascade vertex
 541
 542     esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);
 543     esdCascade->GetPosCovXi(covVtx);
 544     chi2 = esdCascade->GetChi2Xi();
 545
 546     AliAODVertex *vCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex( pos,
 547                                                                      covVtx,
 548                                                                      chi2, // FIXME = Chi2/NDF will be needed
 549                                                                      fPrimaryVertex,
 550                                                                      nCascade, // id
 551                                                                      AliAODVertex::kCascade);
 552     fPrimaryVertex->AddDaughter(vCascade);
 553
 554 //    if (fDebug > 2) {
 555 //      printf("---- Cascade / Cascade Vertex (AOD) : \n");
 556 //      vCascade->Print();
 557 //    }
 558
 559 //    if(esd.GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(const_cast<AliESDEvent*>(&esd), (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production starting form LHC10e without Tender.
 560
 561
 562     // 3 - Add the bachelor track from the cascade
 563
 564     if (!fUsedTrack[idxBachFromCascade]) {
 565
 566       esdCascadeBach->GetPxPyPz(momBach);
 567       esdCascadeBach->GetXYZ(pos);
 568       esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);
 569       esdCascadeBach->GetESDpid(pid);
 570
 571             fUsedTrack[idxBachFromCascade] = kTRUE;
 572             UInt_t selectInfo = 0;
 573             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);
 574             if (fMChandler) fMChandler->SelectParticle(esdCascadeBach->GetLabel());
 575             aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(esdCascadeBach->GetID(),
 576                                                                            esdCascadeBach->GetLabel(),
 577                                                                            momBach,
 578                                                                            kTRUE,
 579                                                                            pos,
 580                                                                            kFALSE, // Why kFALSE for "isDCA" ? FIXME
 581                                                                            covTr,
 582                                                                            (Short_t)esdCascadeBach->GetSign(),
 583                                                                            esdCascadeBach->GetITSClusterMap(),
 584                                                                            pid,
 585                                                                            vCascade,
 586                                                                            kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 587                                                                            vtx->UsesTrack(esdCascadeBach->GetID()),
 588                                                                            AliAODTrack::kSecondary,
 589                                                                            selectInfo);
 590             aodTrack->SetTPCFitMap(esdCascadeBach->GetTPCFitMap());
 591             aodTrack->SetTPCClusterMap(esdCascadeBach->GetTPCClusterMap());
 592             aodTrack->SetTPCSharedMap (esdCascadeBach->GetTPCSharedMap());
 593             aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeBach));
 594             aodTrack->SetTPCPointsF(esdCascadeBach->GetTPCNclsF());
 595             aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeBach->GetTPCCrossedRows()));
 596             fAODTrackRefs->AddAt(aodTrack,idxBachFromCascade);
 597
 598             if (esdCascadeBach->GetSign() > 0) ++fNumberOfPositiveTracks;
 599             aodTrack->ConvertAliPIDtoAODPID();
 600             aodTrack->SetFlags(esdCascadeBach->GetStatus());
 601       SetAODPID(esdCascadeBach,aodTrack,detpid);
 602     }
 603     else {
 604             aodTrack = static_cast<AliAODTrack*>( fAODTrackRefs->At(idxBachFromCascade) );
 605     }
 606
 607     vCascade->AddDaughter(aodTrack);
 608
 609 //    if (fDebug > 4) {
 610 //      printf("---- Cascade / bach dghter : \n");
 611 //      aodTrack->Print();
 612 //    }
 613
 614
 615     // 4 - Add the V0 from the cascade.
 616     // = V0vtx + both pos and neg daughter tracks + the aodV0 itself
 617     //
 618
 619     if ( !fUsedV0[idxV0FromCascade] ) {
 620       // 4.A - if VO structure hasn't been created yet
 621
 622       // 4.A.1 - Create the V0 vertex of the cascade
 623
 624       esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);
 625       esdV0FromCascade->GetPosCov(covVtx);
 626       chi2 = esdV0FromCascade->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3 ?
 627
 628       vV0FromCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex(pos,
 629                                                                          covVtx,
 630                                                                          chi2,
 631                                                                          vCascade,
 632                                                                          idxV0FromCascade, //id of ESDv0
 633                                                                          AliAODVertex::kV0);
 634       // Note:
 635       //    one V0 can be used by several cascades.
 636       // So, one AOD V0 vtx can have several parent vtx.
 637       // This is not directly allowed by AliAODvertex.
 638       // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash
 639       // but to a problem of consistency within AODEvent.
 640       // -> See below paragraph 4.B, for the proposed treatment of such a case.
 641
 642       // Add the vV0FromCascade to the aodVOVtxRefs
 643       fAODV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);
 644
 645
 646       // 4.A.2 - Add the positive tracks from the V0
 647
 648       esdCascadePos->GetPxPyPz(momPos);
 649       esdCascadePos->GetXYZ(pos);
 650       esdCascadePos->GetCovarianceXYZPxPyPz(covTr);
 651       esdCascadePos->GetESDpid(pid);
 652
 653
 654       if (!fUsedTrack[idxPosFromV0Dghter]) {
 655         fUsedTrack[idxPosFromV0Dghter] = kTRUE;
 656
 657         UInt_t selectInfo = 0;
 658         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);
 659         if(fMChandler) fMChandler->SelectParticle(esdCascadePos->GetLabel());
 660         aodTrack = new(tracksArray[fNumberOfTracks++])
 661         AliAODTrack(  esdCascadePos->GetID(),
 662                     esdCascadePos->GetLabel(),
 663                     momPos,
 664                     kTRUE,
 665                     pos,
 666                     kFALSE, // Why kFALSE for "isDCA" ? FIXME
 667                     covTr,
 668                     (Short_t)esdCascadePos->GetSign(),
 669                     esdCascadePos->GetITSClusterMap(),
 670                     pid,
 671                     vV0FromCascade,
 672                     kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 673                     vtx->UsesTrack(esdCascadePos->GetID()),
 674                     AliAODTrack::kSecondary,
 675                     selectInfo);
 676         aodTrack->SetTPCFitMap(esdCascadePos->GetTPCFitMap());
 677         aodTrack->SetTPCClusterMap(esdCascadePos->GetTPCClusterMap());
 678         aodTrack->SetTPCSharedMap (esdCascadePos->GetTPCSharedMap());
 679         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadePos));
 680         aodTrack->SetTPCPointsF(esdCascadePos->GetTPCNclsF());
 681         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadePos->GetTPCCrossedRows()));
 682         fAODTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);
 683
 684         if (esdCascadePos->GetSign() > 0) ++fNumberOfPositiveTracks;
 685         aodTrack->ConvertAliPIDtoAODPID();
 686         aodTrack->SetFlags(esdCascadePos->GetStatus());
 687         SetAODPID(esdCascadePos,aodTrack,detpid);
 688       }
 689       else {
 690         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxPosFromV0Dghter));
 691       }
 692       vV0FromCascade->AddDaughter(aodTrack);
 693
 694
 695       // 4.A.3 - Add the negative tracks from the V0
 696
 697       esdCascadeNeg->GetPxPyPz(momNeg);
 698       esdCascadeNeg->GetXYZ(pos);
 699       esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);
 700       esdCascadeNeg->GetESDpid(pid);
 701
 702
 703       if (!fUsedTrack[idxNegFromV0Dghter]) {
 704         fUsedTrack[idxNegFromV0Dghter] = kTRUE;
 705
 706         UInt_t selectInfo = 0;
 707         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);
 708         if(fMChandler)fMChandler->SelectParticle(esdCascadeNeg->GetLabel());
 709         aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(  esdCascadeNeg->GetID(),
 710                                                       esdCascadeNeg->GetLabel(),
 711                                                       momNeg,
 712                                                       kTRUE,
 713                                                       pos,
 714                                                       kFALSE, // Why kFALSE for "isDCA" ? FIXME
 715                                                       covTr,
 716                                                       (Short_t)esdCascadeNeg->GetSign(),
 717                                                       esdCascadeNeg->GetITSClusterMap(),
 718                                                       pid,
 719                                                       vV0FromCascade,
 720                                                       kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 721                                                       vtx->UsesTrack(esdCascadeNeg->GetID()),
 722                                                       AliAODTrack::kSecondary,
 723                                                       selectInfo);
 724         aodTrack->SetTPCFitMap(esdCascadeNeg->GetTPCFitMap());
 725         aodTrack->SetTPCClusterMap(esdCascadeNeg->GetTPCClusterMap());
 726         aodTrack->SetTPCSharedMap (esdCascadeNeg->GetTPCSharedMap());
 727         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeNeg));
 728         aodTrack->SetTPCPointsF(esdCascadeNeg->GetTPCNclsF());
 729         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeNeg->GetTPCCrossedRows()));
 730         fAODTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);
 731
 732         if (esdCascadeNeg->GetSign() > 0) ++fNumberOfPositiveTracks;
 733         aodTrack->ConvertAliPIDtoAODPID();
 734         aodTrack->SetFlags(esdCascadeNeg->GetStatus());
 735         SetAODPID(esdCascadeNeg,aodTrack,detpid);
 736       }
 737       else {
 738         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxNegFromV0Dghter));
 739       }
 740
 741       vV0FromCascade->AddDaughter(aodTrack);
 742
 743
 744       // 4.A.4 - Add the V0 from cascade to the V0 array
 745
 746       Double_t  dcaV0Daughters      = esdV0FromCascade->GetDcaV0Daughters();
 747       Double_t  dcaV0ToPrimVertex   = esdV0FromCascade->GetD( esd.GetPrimaryVertex()->GetX(),
 748                                                              esd.GetPrimaryVertex()->GetY(),
 749                                                              esd.GetPrimaryVertex()->GetZ() );
 750       esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] );
 751       esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] );
 752
 753       Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
 754       dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD(      esd.GetPrimaryVertex()->GetX(),
 755                                                                   esd.GetPrimaryVertex()->GetY(),
 756                                                                   esd.GetMagneticField())        );
 757       dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD(      esd.GetPrimaryVertex()->GetX(),
 758                                                                   esd.GetPrimaryVertex()->GetY(),
 759                                                                   esd.GetMagneticField())        );
 760
 761       aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0( vV0FromCascade,
 762                                                   dcaV0Daughters,
 763                                                   dcaV0ToPrimVertex,
 764                                                   momPosAtV0vtx,
 765                                                   momNegAtV0vtx,
 766                                                   dcaDaughterToPrimVertex);
 767       // set the aod v0 on-the-fly status
 768       aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());
 769
 770       // Add the aodV0 to the aodVORefs
 771       fAODV0Refs->AddAt(aodV0,idxV0FromCascade);
 772
 773       fUsedV0[idxV0FromCascade] = kTRUE;
 774
 775     } else {
 776       // 4.B - if V0 structure already used
 777
 778       // Note :
 779       //    one V0 can be used by several cascades (frequent in PbPb evts) :
 780       // same V0 which used but attached to different bachelor tracks
 781       // -> aodVORefs and fAODV0VtxRefs are needed.
 782       // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.
 783
 784       vV0FromCascade = static_cast<AliAODVertex*>( fAODV0VtxRefs->At(idxV0FromCascade) );
 785       aodV0          = static_cast<AliAODv0*>    ( fAODV0Refs   ->At(idxV0FromCascade) );
 786
 787       // - Treatment of the parent for such a "re-used" V0 :
 788       // Insert the cascade that reuses the V0 vertex in the lineage chain
 789       // Before : vV0 -> vCascade1 -> vPrimary
 790       //  - Hyp : cascade2 uses the same V0 as cascade1
 791       //  After :  vV0 -> vCascade2 -> vCascade1 -> vPrimary
 792
 793       AliAODVertex *vCascadePreviousParent = static_cast<AliAODVertex*> (vV0FromCascade->GetParent());
 794                         vV0FromCascade->SetParent(vCascade);
 795                         vCascade      ->SetParent(vCascadePreviousParent);
 796
 797 //      if(fDebug > 2)
 798 //        printf("---- Cascade / Lineage insertion\n"
 799 //               "Parent of V0 vtx                 = Cascade vtx %p\n"
 800 //               "Parent of the cascade vtx        = Cascade vtx %p\n"
 801 //               "Parent of the parent cascade vtx = Cascade vtx %p\n",
 802 //               static_cast<void*> (vV0FromCascade->GetParent()),
 803 //               static_cast<void*> (vCascade->GetParent()),
 804 //               static_cast<void*> (vCascadePreviousParent->GetParent()) );
 805
 806     }// end if V0 structure already used
 807
 808 //    if (fDebug > 2) {
 809 //      printf("---- Cascade / V0 vertex: \n");
 810 //      vV0FromCascade->Print();
 811 //    }
 812 //
 813 //    if (fDebug > 4) {
 814 //      printf("---- Cascade / pos dghter : \n");
 815 //                      aodTrack->Print();
 816 //      printf("---- Cascade / neg dghter : \n");
 817 //                      aodTrack->Print();
 818 //      printf("---- Cascade / aodV0 : \n");
 819 //                      aodV0->Print();
 820 //    }
 821
 822     // In any case (used V0 or not), add the V0 vertex to the cascade one.
 823     vCascade->AddDaughter(vV0FromCascade);
 824
 825
 826     // 5 - Add the primary track of the cascade (if any)
 827
 828
 829     // 6 - Add the cascade to the AOD array of cascades
 830
 831     Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd.GetPrimaryVertex()->GetX(),
 832                                                                      esd.GetPrimaryVertex()->GetY(),
 833                                                                      esd.GetMagneticField())        );
 834
 835     Double_t momBachAtCascadeVtx[3]={0.};
 836
 837     esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);
 838
 839     aodCascade = new(cascadesArray[fNumberOfCascades++]) AliAODcascade(  vCascade,
 840                                                           esdCascade->Charge(),
 841                                                           esdCascade->GetDcaXiDaughters(),
 842                                                           -999.,
 843                                                           // DCAXiToPrimVtx -> needs to be calculated   ----|
 844                                                           // doesn't exist at ESD level;
 845                                                           // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)
 846                                                           dcaBachToPrimVertexXY,
 847                                                           momBachAtCascadeVtx,
 848                                                           *aodV0);
 849
 850     if (fDebug > 10) {
 851       printf("---- Cascade / AOD cascade : \n\n");
 852       aodCascade->PrintXi(fPrimaryVertex->GetX(), fPrimaryVertex->GetY(), fPrimaryVertex->GetZ());
 853     }
 854
 855   } // end of the loop on cascades
 856
 857   Cascades().Expand(fNumberOfCascades);
 858 }
 859
 860 //______________________________________________________________________________
 861 void AliAnalysisTaskESDfilter::ConvertV0s(const AliESDEvent& esd)
 862 {
 863   // Access to the AOD container of V0s
 864
 865   AliCodeTimerAuto("",0);
 866
 867   //
 868   // V0s
 869   //
 870
 871   Double_t pos[3] = { 0. };
 872   Double_t chi2(0.0);
 873   Double_t covVtx[6] = { 0. };
 874   Double_t momPos[3]={0.};
 875   Double_t covTr[21]={0.};
 876   Double_t pid[10]={0.};
 877   AliAODTrack* aodTrack(0x0);
 878   AliAODPid* detpid(0x0);
 879   Double_t momNeg[3]={0.};
 880   Double_t momPosAtV0vtx[3]={0.};
 881   Double_t momNegAtV0vtx[3]={0.};
 882
 883   for (Int_t nV0 = 0; nV0 < esd.GetNumberOfV0s(); ++nV0)
 884   {
 885     if (fUsedV0[nV0]) continue; // skip if already added to the AOD
 886
 887     AliESDv0 *v0 = esd.GetV0(nV0);
 888     Int_t posFromV0 = v0->GetPindex();
 889     Int_t negFromV0 = v0->GetNindex();
 890
 891     // V0 selection
 892     //
 893     AliESDVertex *esdVtx   = new AliESDVertex(*(esd.GetPrimaryVertex()));
 894     AliESDtrack  *esdV0Pos = esd.GetTrack(posFromV0);
 895     AliESDtrack  *esdV0Neg = esd.GetTrack(negFromV0);
 896     TList v0objects;
 897     v0objects.AddAt(v0,                      0);
 898     v0objects.AddAt(esdV0Pos,                1);
 899     v0objects.AddAt(esdV0Neg,                2);
 900     v0objects.AddAt(esdVtx,                  3);
 901     UInt_t selectV0 = 0;
 902     if (fV0Filter) {
 903       selectV0 = fV0Filter->IsSelected(&v0objects);
 904       // this is a little awkward but otherwise the
 905       // list wants to access the pointer (delete it)
 906       // again when going out of scope
 907       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 908       esdVtx = 0;
 909       if (!selectV0)
 910         continue;
 911     }
 912     else{
 913       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 914       esdVtx = 0;
 915     }
 916
 917     v0->GetXYZ(pos[0], pos[1], pos[2]);
 918
 919     if (!fOldESDformat) {
 920             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
 921             v0->GetPosCov(covVtx);
 922     } else {
 923             chi2 = -999.;
 924             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 925     }
 926
 927
 928     AliAODVertex * vV0 =
 929           new(Vertices()[fNumberOfVertices++]) AliAODVertex(pos,
 930                                                       covVtx,
 931                                                       chi2,
 932                                                       fPrimaryVertex,
 933                                                       nV0,
 934                                                       AliAODVertex::kV0);
 935     fPrimaryVertex->AddDaughter(vV0);
 936
 937
 938     // Add the positive tracks from the V0
 939
 940
 941     esdV0Pos->GetPxPyPz(momPos);
 942     esdV0Pos->GetXYZ(pos);
 943     esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
 944     esdV0Pos->GetESDpid(pid);
 945
 946     const AliESDVertex *vtx = esd.GetPrimaryVertex();
 947
 948     if (!fUsedTrack[posFromV0]) {
 949             fUsedTrack[posFromV0] = kTRUE;
 950             UInt_t selectInfo = 0;
 951             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);
 952             if(fMChandler)fMChandler->SelectParticle(esdV0Pos->GetLabel());
 953             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Pos->GetID(),
 954                                                     esdV0Pos->GetLabel(),
 955                                                     momPos,
 956                                                     kTRUE,
 957                                                     pos,
 958                                                     kFALSE,
 959                                                     covTr,
 960                                                     (Short_t)esdV0Pos->GetSign(),
 961                                                     esdV0Pos->GetITSClusterMap(),
 962                                                     pid,
 963                                                     vV0,
 964                                                     kTRUE,  // check if this is right
 965                                                     vtx->UsesTrack(esdV0Pos->GetID()),
 966                                                     AliAODTrack::kSecondary,
 967                                                     selectInfo);
 968             aodTrack->SetTPCFitMap(esdV0Pos->GetTPCFitMap());
 969             aodTrack->SetTPCClusterMap(esdV0Pos->GetTPCClusterMap());
 970             aodTrack->SetTPCSharedMap (esdV0Pos->GetTPCSharedMap());
 971             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Pos));
 972             aodTrack->SetTPCPointsF(esdV0Pos->GetTPCNclsF());
 973             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Pos->GetTPCCrossedRows()));
 974             fAODTrackRefs->AddAt(aodTrack,posFromV0);
 975             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
 976             if (esdV0Pos->GetSign() > 0) ++fNumberOfPositiveTracks;
 977             aodTrack->ConvertAliPIDtoAODPID();
 978             aodTrack->SetFlags(esdV0Pos->GetStatus());
 979       SetAODPID(esdV0Pos,aodTrack,detpid);
 980     }
 981     else {
 982             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));
 983             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
 984     }
 985     vV0->AddDaughter(aodTrack);
 986
 987     // Add the negative tracks from the V0
 988
 989     esdV0Neg->GetPxPyPz(momNeg);
 990     esdV0Neg->GetXYZ(pos);
 991     esdV0Neg->GetCovarianceXYZPxPyPz(covTr);
 992     esdV0Neg->GetESDpid(pid);
 993
 994     if (!fUsedTrack[negFromV0]) {
 995             fUsedTrack[negFromV0] = kTRUE;
 996             UInt_t selectInfo = 0;
 997             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);
 998             if(fMChandler)fMChandler->SelectParticle(esdV0Neg->GetLabel());
 999             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Neg->GetID(),
1000                                                     esdV0Neg->GetLabel(),
1001                                                     momNeg,
1002                                                     kTRUE,
1003                                                     pos,
1004                                                     kFALSE,
1005                                                     covTr,
1006                                                     (Short_t)esdV0Neg->GetSign(),
1007                                                     esdV0Neg->GetITSClusterMap(),
1008                                                     pid,
1009                                                     vV0,
1010                                                     kTRUE,  // check if this is right
1011                                                     vtx->UsesTrack(esdV0Neg->GetID()),
1012                                                     AliAODTrack::kSecondary,
1013                                                     selectInfo);
1014             aodTrack->SetTPCFitMap(esdV0Neg->GetTPCFitMap());
1015             aodTrack->SetTPCClusterMap(esdV0Neg->GetTPCClusterMap());
1016             aodTrack->SetTPCSharedMap (esdV0Neg->GetTPCSharedMap());
1017             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Neg));
1018             aodTrack->SetTPCPointsF(esdV0Neg->GetTPCNclsF());
1019             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Neg->GetTPCCrossedRows()));
1020
1021             fAODTrackRefs->AddAt(aodTrack,negFromV0);
1022             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
1023             if (esdV0Neg->GetSign() > 0) ++fNumberOfPositiveTracks;
1024             aodTrack->ConvertAliPIDtoAODPID();
1025             aodTrack->SetFlags(esdV0Neg->GetStatus());
1026       SetAODPID(esdV0Neg,aodTrack,detpid);
1027     }
1028     else {
1029             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(negFromV0));
1030             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
1031     }
1032     vV0->AddDaughter(aodTrack);
1033
1034
1035     // Add the V0 the V0 array as well
1036
1037     Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();
1038     Double_t  dcaV0ToPrimVertex   = v0->GetD(esd.GetPrimaryVertex()->GetX(),
1039                                              esd.GetPrimaryVertex()->GetY(),
1040                                              esd.GetPrimaryVertex()->GetZ());
1041
1042     v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]);
1043     v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]);
1044
1045     Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
1046     dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD(  esd.GetPrimaryVertex()->GetX(),
1047                                                            esd.GetPrimaryVertex()->GetY(),
1048                                                            esd.GetMagneticField()) );
1049     dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD(  esd.GetPrimaryVertex()->GetX(),
1050                                                            esd.GetPrimaryVertex()->GetY(),
1051                                                            esd.GetMagneticField()) );
1052
1053     AliAODv0* aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0(vV0,
1054                                                 dcaV0Daughters,
1055                                                 dcaV0ToPrimVertex,
1056                                                 momPosAtV0vtx,
1057                                                 momNegAtV0vtx,
1058                                                 dcaDaughterToPrimVertex);
1059
1060     // set the aod v0 on-the-fly status
1061     aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
1062   }//End of loop on V0s
1063
1064   V0s().Expand(fNumberOfV0s);
1065 }
1066
1067 //______________________________________________________________________________
1068 void AliAnalysisTaskESDfilter::ConvertTPCOnlyTracks(const AliESDEvent& esd)
1069 {
1070   // Convert TPC only tracks
1071   // Here we have wo hybrid appraoch to remove fakes
1072   // ******* ITSTPC ********
1073   // Uses a cut on the ITS properties to select global tracks
1074   // which are than marked as HybdridITSTPC for the remainder
1075   // the TPC only tracks are flagged as HybridITSTPConly.
1076   // Note, in order not to get fakes back in the TPC cuts, one needs
1077   // two "ITS" cuts one tight (1) (to throw out fakes) and one lose (2) (to NOT flag the trakcs in the TPC only)
1078   // using cut number (3)
1079   // so fHybridFilterMask == (1)|(2) fTPCFilterMask = (3), Usercode needs to slect with mask = (1)|(3) and track->IsHybridITSTPC()
1080   // ******* TPC ********
1081   // Here, only TPC tracks are flagged that pass the tight ITS cuts and tracks that pass the TPC cuts and NOT the loose ITS cuts
1082   // the ITS cuts neeed to be added to the filter as extra cuts, since here the selections info is reset in the global and put to the TPC only track
1083
1084   AliCodeTimerAuto("",0);
1085
1086   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1087   for(int it = 0;it < fNumberOfTracks;++it)
1088   {
1089     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1090     if(!tr)continue;
1091     UInt_t map = tr->GetFilterMap();
1092     if(map&fTPCConstrainedFilterMask){
1093       // we only reset the track select ionfo, no deletion...
1094       tr->SetFilterMap(map&~fTPCConstrainedFilterMask);
1095     }
1096     if(map&fHybridFilterMaskTPCCG){
1097       // this is one part of the hybrid tracks
1098       // the others not passing the selection will be TPC only selected below
1099       tr->SetIsHybridTPCConstrainedGlobal(kTRUE);
1100     }
1101   }
1102   // Loop over the ESD trcks and pick out the tracks passing TPC only cuts
1103
1104
1105   const AliESDVertex *vtxSPD = esd.GetPrimaryVertexSPD();
1106   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1107
1108   Double_t pos[3] = { 0. };
1109   Double_t covTr[21]={0.};
1110   Double_t pid[10]={0.};
1111
1112
1113   Double_t p[3] = { 0. };
1114
1115   Double_t pDCA[3] = { 0. }; // momentum at DCA
1116   Double_t rDCA[3] = { 0. }; // position at DCA
1117   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1118   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1119
1120
1121   AliAODTrack* aodTrack(0x0);
1122   //  AliAODPid* detpid(0x0);
1123
1124   // account for change in pT after the constraint
1125   Float_t ptMax = 1E10;
1126   Float_t ptMin = 0;
1127   for(int i = 0;i<32;i++){
1128     if(fTPCConstrainedFilterMask&(1<<i)){
1129       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1130       Float_t tmp1= 0,tmp2 = 0;
1131       cuts->GetPtRange(tmp1,tmp2);
1132       if(tmp1>ptMin)ptMin=tmp1;
1133       if(tmp2<ptMax)ptMax=tmp2;
1134     }
1135   }
1136
1137   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1138   {
1139     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1140
1141     UInt_t selectInfo = 0;
1142     Bool_t isHybridITSTPC = false;
1143     //
1144     // Track selection
1145     if (fTrackFilter) {
1146       selectInfo = fTrackFilter->IsSelected(esdTrack);
1147     }
1148
1149     if(!(selectInfo&fHybridFilterMaskTPCCG)){
1150       // not already selected tracks, use second part of hybrid tracks
1151       isHybridITSTPC = true;
1152       // too save space one could only store these...
1153     }
1154
1155     selectInfo &= fTPCConstrainedFilterMask;
1156     if (!selectInfo)continue;
1157     if (fWriteHybridTPCCOnly&&!isHybridITSTPC)continue; // write only complementary tracks
1158     // create a tpc only tracl
1159     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(&esd),esdTrack->GetID());
1160     if(!track) continue;
1161
1162     if(track->Pt()>0.)
1163     {
1164       // only constrain tracks above threshold
1165       AliExternalTrackParam exParam;
1166       // take the B-field from the ESD, no 3D fieldMap available at this point
1167       Bool_t relate = false;
1168       relate = track->RelateToVertexTPC(vtxSPD,esd.GetMagneticField(),kVeryBig,&exParam);
1169       if(!relate){
1170         delete track;
1171         continue;
1172       }
1173       // fetch the track parameters at the DCA (unconstraint)
1174       if(track->GetTPCInnerParam()){
1175         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1176         track->GetTPCInnerParam()->GetXYZ(rDCA);
1177       }
1178       // get the DCA to the vertex:
1179       track->GetImpactParametersTPC(dDCA,cDCA);
1180       // set the constrained parameters to the track
1181       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1182     }
1183
1184     track->GetPxPyPz(p);
1185
1186     Float_t pT = track->Pt();
1187     if(pT<ptMin||pT>ptMax){
1188       delete track;
1189       continue;
1190     }
1191
1192     //
1193
1194
1195     track->GetXYZ(pos);
1196     track->GetCovarianceXYZPxPyPz(covTr);
1197     esdTrack->GetESDpid(pid);// original PID
1198
1199     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1200     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((track->GetID()+1)*-1,
1201                                                             track->GetLabel(),
1202                                                             p,
1203                                                             kTRUE,
1204                                                             pos,
1205                                                             kFALSE,
1206                                                             covTr,
1207                                                             (Short_t)track->GetSign(),
1208                                                             track->GetITSClusterMap(),
1209                                                             pid,
1210                                                             fPrimaryVertex,
1211                                                             kTRUE, // check if this is right
1212                                                             vtx->UsesTrack(track->GetID()),
1213                                                             AliAODTrack::kPrimary,
1214                                                             selectInfo);
1215     aodTrack->SetIsHybridTPCConstrainedGlobal(isHybridITSTPC);
1216     aodTrack->SetTPCFitMap(track->GetTPCFitMap());
1217     aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
1218     aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
1219     aodTrack->SetIsTPCConstrained(kTRUE);
1220     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack)); // original track
1221     // set the DCA values to the AOD track
1222     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1223     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1224     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1225
1226     aodTrack->SetFlags(track->GetStatus());
1227     aodTrack->SetTPCPointsF(track->GetTPCNclsF());
1228     aodTrack->SetTPCNCrossedRows(UShort_t(track->GetTPCCrossedRows()));
1229
1230     //Perform progagation of tracks if needed
1231     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1232     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1233
1234     // do not duplicate PID information
1235     //    aodTrack->ConvertAliPIDtoAODPID();
1236     //    SetAODPID(esdTrack,aodTrack,detpid);
1237
1238     delete track;
1239   } // end of loop on tracks
1240
1241 }
1242
1243 //______________________________________________________________________________
1244 void AliAnalysisTaskESDfilter::ConvertGlobalConstrainedTracks(const AliESDEvent& esd)
1245 {
1246
1247   // Here we have the option to store the complement from global constraint information
1248   // to tracks passing tight cuts (1) in order not to get fakes back in, one needs
1249   // two sets of cuts one tight (1) (to throw out fakes) and one lose (2) (fakes/bad tracks would pass (2) but not (1))
1250   // using cut number (3) selects the tracks that complement (1) e.g. tracks witout ITS refit or cluster requirement
1251
1252
1253   AliCodeTimerAuto("",0);
1254
1255   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1256   for(int it = 0;it < fNumberOfTracks;++it)
1257   {
1258     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1259     if(!tr)continue;
1260     UInt_t map = tr->GetFilterMap();
1261     if(map&fGlobalConstrainedFilterMask){
1262       // we only reset the track select info, no deletion...
1263       // mask reset mask in case track is already taken
1264       tr->SetFilterMap(map&~fGlobalConstrainedFilterMask);
1265     }
1266     if(map&fHybridFilterMaskGCG){
1267       // this is one part of the hybrid tracks
1268       // the others not passing the selection will be the ones selected below
1269       tr->SetIsHybridGlobalConstrainedGlobal(kTRUE);
1270     }
1271   }
1272   // Loop over the ESD trcks and pick out the tracks passing the GlobalConstraint cuts
1273
1274
1275   Double_t pos[3] = { 0. };
1276   Double_t covTr[21]={0.};
1277   Double_t pid[10]={0.};
1278   Double_t p[3] = { 0. };
1279
1280   Double_t pDCA[3] = { 0. }; // momentum at DCA
1281   Double_t rDCA[3] = { 0. }; // position at DCA
1282   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1283   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1284
1285
1286   AliAODTrack* aodTrack(0x0);
1287   AliAODPid* detpid(0x0);
1288   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1289
1290   // account for change in pT after the constraint
1291   Float_t ptMax = 1E10;
1292   Float_t ptMin = 0;
1293   for(int i = 0;i<32;i++){
1294     if(fGlobalConstrainedFilterMask&(1<<i)){
1295       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1296       Float_t tmp1= 0,tmp2 = 0;
1297       cuts->GetPtRange(tmp1,tmp2);
1298       if(tmp1>ptMin)ptMin=tmp1;
1299       if(tmp2<ptMax)ptMax=tmp2;
1300     }
1301   }
1302
1303
1304
1305  for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1306   {
1307     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1308     const AliExternalTrackParam * exParamGC = esdTrack->GetConstrainedParam();
1309     if(!exParamGC)continue;
1310
1311     UInt_t selectInfo = 0;
1312     Bool_t isHybridGC = false;
1313
1314     //
1315     // Track selection
1316     if (fTrackFilter) {
1317       selectInfo = fTrackFilter->IsSelected(esdTrack);
1318     }
1319
1320
1321     if(!(selectInfo&fHybridFilterMaskGCG))isHybridGC = true;
1322     if (fWriteHybridGCOnly&&!isHybridGC)continue; // write only complementary tracks
1323
1324     selectInfo &= fGlobalConstrainedFilterMask;
1325     if (!selectInfo)continue;
1326     // fetch the track parameters at the DCA (unconstrained)
1327     esdTrack->GetPxPyPz(pDCA);
1328     esdTrack->GetXYZ(rDCA);
1329     // get the DCA to the vertex:
1330     esdTrack->GetImpactParameters(dDCA,cDCA);
1331
1332     if (!esdTrack->GetConstrainedPxPyPz(p)) continue;
1333
1334
1335     Float_t pT = exParamGC->Pt();
1336     if(pT<ptMin||pT>ptMax){
1337       continue;
1338     }
1339
1340
1341     esdTrack->GetConstrainedXYZ(pos);
1342     exParamGC->GetCovarianceXYZPxPyPz(covTr);
1343     esdTrack->GetESDpid(pid);
1344     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1345     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((esdTrack->GetID()+1)*-1,
1346                                                             esdTrack->GetLabel(),
1347                                                             p,
1348                                                             kTRUE,
1349                                                             pos,
1350                                                             kFALSE,
1351                                                             covTr,
1352                                                             (Short_t)esdTrack->GetSign(),
1353                                                             esdTrack->GetITSClusterMap(),
1354                                                             pid,
1355                                                             fPrimaryVertex,
1356                                                             kTRUE, // check if this is right
1357                                                             vtx->UsesTrack(esdTrack->GetID()),
1358                                                             AliAODTrack::kPrimary,
1359                                                             selectInfo);
1360     aodTrack->SetIsHybridGlobalConstrainedGlobal(isHybridGC);
1361     aodTrack->SetIsGlobalConstrained(kTRUE);
1362     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1363     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1364     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1365     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1366
1367
1368     // set the DCA values to the AOD track
1369     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1370     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1371     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1372
1373     aodTrack->SetFlags(esdTrack->GetStatus());
1374     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1375     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1376
1377     if(isHybridGC){
1378       // only copy AOD information for hybrid, no duplicate information
1379       aodTrack->ConvertAliPIDtoAODPID();
1380       SetAODPID(esdTrack,aodTrack,detpid);
1381     }
1382
1383     //Perform progagation of tracks if needed
1384     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1385     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1386   } // end of loop on tracks
1387
1388 }
1389
1390
1391 //______________________________________________________________________________
1392 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)
1393 {
1394   // Tracks (primary and orphan)
1395
1396   AliCodeTimerAuto("",0);
1397
1398   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));
1399
1400   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1401   Double_t p[3] = { 0. };
1402   Double_t pos[3] = { 0. };
1403   Double_t covTr[21] = { 0. };
1404   Double_t pid[10] = { 0. };
1405   AliAODTrack* aodTrack(0x0);
1406   AliAODPid* detpid(0x0);
1407
1408   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1409   {
1410     if (fUsedTrack[nTrack]) continue;
1411
1412     AliESDtrack *esdTrack = esd.GetTrack(nTrack);
1413     UInt_t selectInfo = 0;
1414     //
1415     // Track selection
1416     if (fTrackFilter) {
1417             selectInfo = fTrackFilter->IsSelected(esdTrack);
1418             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;
1419     }
1420
1421
1422     esdTrack->GetPxPyPz(p);
1423     esdTrack->GetXYZ(pos);
1424     esdTrack->GetCovarianceXYZPxPyPz(covTr);
1425     esdTrack->GetESDpid(pid);
1426     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1427     fPrimaryVertex->AddDaughter(aodTrack =
1428                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),
1429                                                             esdTrack->GetLabel(),
1430                                                             p,
1431                                                             kTRUE,
1432                                                             pos,
1433                                                             kFALSE,
1434                                                             covTr,
1435                                                             (Short_t)esdTrack->GetSign(),
1436                                                             esdTrack->GetITSClusterMap(),
1437                                                             pid,
1438                                                             fPrimaryVertex,
1439                                                             kTRUE, // check if this is right
1440                                                             vtx->UsesTrack(esdTrack->GetID()),
1441                                                             AliAODTrack::kPrimary,
1442                                                             selectInfo)
1443                          );
1444     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1445     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1446     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1447     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1448     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1449     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1450     if(esdTrack->IsEMCAL()) aodTrack->SetEMCALcluster(esdTrack->GetEMCALcluster());
1451     if(esdTrack->IsPHOS())  aodTrack->SetPHOScluster(esdTrack->GetPHOScluster());
1452
1453     //Perform progagation of tracks if needed
1454     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1455     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1456
1457     fAODTrackRefs->AddAt(aodTrack, nTrack);
1458
1459
1460     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;
1461     aodTrack->SetFlags(esdTrack->GetStatus());
1462     aodTrack->ConvertAliPIDtoAODPID();
1463     SetAODPID(esdTrack,aodTrack,detpid);
1464   } // end of loop on tracks
1465 }
1466
1467 //______________________________________________________________________________
1468 void AliAnalysisTaskESDfilter::PropagateTrackToEMCal(AliESDtrack *esdTrack)
1469 {
1470   Double_t trkPos[3] = {0.,0.,0.};
1471   Double_t EMCalEta=-999, EMCalPhi=-999;
1472   Double_t trkphi = esdTrack->Phi()*TMath::RadToDeg();
1473   if(TMath::Abs(esdTrack->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
1474     {
1475       AliExternalTrackParam *trkParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1476       if(trkParam)
1477         {
1478           AliExternalTrackParam trkParamTmp(*trkParam);
1479           if(AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, 430, esdTrack->GetMass(), 20, kTRUE, 0.8, -1))
1480             {
1481               trkParamTmp.GetXYZ(trkPos);
1482               TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1483               EMCalEta = trkPosVec.Eta();
1484               EMCalPhi = trkPosVec.Phi();
1485               if(EMCalPhi<0)  EMCalPhi += 2*TMath::Pi();
1486               esdTrack->SetTrackPhiEtaOnEMCal(EMCalPhi,EMCalEta);
1487             }
1488         }
1489     }
1490 }
1491
1492 //______________________________________________________________________________
1493 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)
1494 {
1495 // Convert PMD Clusters
1496   AliCodeTimerAuto("",0);
1497   Int_t jPmdClusters=0;
1498   // Access to the AOD container of PMD clusters
1499   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
1500   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {
1501     // file pmd clusters, to be revised!
1502     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);
1503     Int_t nLabel = 0;
1504     Int_t *label = 0x0;
1505     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
1506     Double_t pidPmd[13] = { 0.}; // to be revised!
1507     // type not set!
1508     // assoc cluster not set
1509     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
1510   }
1511 }
1512
1513
1514 //______________________________________________________________________________
1515 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)
1516 {
1517 // Convert Calorimeter Clusters
1518   AliCodeTimerAuto("",0);
1519
1520   // Access to the AOD container of clusters
1521   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
1522   Int_t jClusters(0);
1523
1524   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
1525
1526     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
1527
1528     Int_t  id        = cluster->GetID();
1529     Int_t  nLabel    = cluster->GetNLabels();
1530     Int_t *labels    = cluster->GetLabels();
1531     if(labels){
1532                   for(int i = 0;i < nLabel;++i){
1533                           if(fMChandler)fMChandler->SelectParticle(labels[i]);
1534                   }
1535           }
1536
1537     Float_t energy = cluster->E();
1538     Float_t posF[3] = { 0.};
1539     cluster->GetPosition(posF);
1540
1541     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
1542                                                                                       nLabel,
1543                                                                                       labels,
1544                                                                                       energy,
1545                                                                                       posF,
1546                                                                                       NULL,
1547                                                                                       cluster->GetType(),0);
1548
1549     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
1550                                 cluster->GetDispersion(),
1551                                 cluster->GetM20(), cluster->GetM02(),
1552                                 cluster->GetEmcCpvDistance(),
1553                                 cluster->GetNExMax(),cluster->GetTOF()) ;
1554
1555     caloCluster->SetPIDFromESD(cluster->GetPID());
1556     caloCluster->SetNCells(cluster->GetNCells());
1557     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
1558     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
1559
1560     caloCluster->SetTrackDistance(cluster->GetTrackDx(), cluster->GetTrackDz());
1561
1562     Int_t nMatchCount = 0;
1563     TArrayI* matchedT =         cluster->GetTracksMatched();
1564     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
1565       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
1566         Int_t iESDtrack = matchedT->At(im);;
1567         if (fAODTrackRefs->At(iESDtrack) != 0) {
1568           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
1569           nMatchCount++;
1570         }
1571       }
1572     }
1573     if(nMatchCount==0)
1574       caloCluster->SetTrackDistance(-999,-999);
1575
1576   }
1577   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
1578 }
1579
1580 //______________________________________________________________________________
1581 void AliAnalysisTaskESDfilter::ConvertCaloTrigger(TString calo, const AliESDEvent& esd)
1582 {
1583         AliCodeTimerAuto("",0);
1584
1585         if (calo == "PHOS")
1586         {
1587           AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1588           AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1589
1590           aodTrigger.Allocate(esdTrigger.GetEntries());
1591           esdTrigger.Reset();
1592
1593           Float_t a;
1594           Int_t tmod,tabsId;
1595
1596           while (esdTrigger.Next()) {
1597             esdTrigger.GetPosition(tmod,tabsId);
1598             esdTrigger.GetAmplitude(a);
1599             aodTrigger.Add(tmod,tabsId,a,0.,(Int_t*)NULL,0,0,0);
1600           }
1601
1602           return;
1603         }
1604
1605         AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1606
1607         if (aodHandler)
1608         {
1609                 TTree *aodTree = aodHandler->GetTree();
1610
1611                 if (aodTree)
1612                 {
1613                         Int_t *type = esd.GetCaloTriggerType();
1614
1615                         for (Int_t i = 0; i < 15; i++)
1616                         {
1617                                 aodTree->GetUserInfo()->Add(new TParameter<int>(Form("EMCALCaloTrigger%d",i), type[i]));
1618                         }
1619                 }
1620         }
1621
1622         AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1623
1624         AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1625
1626         aodTrigger.Allocate(esdTrigger.GetEntries());
1627
1628         esdTrigger.Reset();
1629         while (esdTrigger.Next())
1630         {
1631                 Int_t px, py, ts, nTimes, times[10], b;
1632                 Float_t a, t;
1633
1634                 esdTrigger.GetPosition(px, py);
1635
1636                 esdTrigger.GetAmplitude(a);
1637                 esdTrigger.GetTime(t);
1638
1639                 esdTrigger.GetL0Times(times);
1640                 esdTrigger.GetNL0Times(nTimes);
1641
1642                 esdTrigger.GetL1TimeSum(ts);
1643
1644                 esdTrigger.GetTriggerBits(b);
1645
1646                 aodTrigger.Add(px, py, a, t, times, nTimes, ts, b);
1647         }
1648
1649         for (int i = 0; i < 4; i++) aodTrigger.SetL1Threshold(i, esdTrigger.GetL1Threshold(i));
1650
1651         Int_t v0[2] =
1652         {
1653                 esdTrigger.GetL1V0(0),
1654                 esdTrigger.GetL1V0(1)
1655         };
1656
1657         aodTrigger.SetL1V0(v0);
1658         aodTrigger.SetL1FrameMask(esdTrigger.GetL1FrameMask());
1659 }
1660
1661 //______________________________________________________________________________
1662 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)
1663 {
1664 // Convert EMCAL Cells
1665   AliCodeTimerAuto("",0);
1666   // fill EMCAL cell info
1667   if (esd.GetEMCALCells()) { // protection against missing ESD information
1668     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
1669     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
1670
1671     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1672     aodEMcells.CreateContainer(nEMcell);
1673     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
1674     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
1675       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),
1676                          esdEMcells.GetTime(iCell), esdEMcells.GetMCLabel(iCell), esdEMcells.GetEFraction(iCell));
1677     }
1678     aodEMcells.Sort();
1679   }
1680 }
1681
1682 //______________________________________________________________________________
1683 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)
1684 {
1685 // Convert PHOS Cells
1686   AliCodeTimerAuto("",0);
1687   // fill PHOS cell info
1688   if (esd.GetPHOSCells()) { // protection against missing ESD information
1689     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
1690     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
1691
1692     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
1693     aodPHcells.CreateContainer(nPHcell);
1694     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
1695     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
1696       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),
1697                          esdPHcells.GetTime(iCell), esdPHcells.GetMCLabel(iCell), esdPHcells.GetEFraction(iCell));
1698     }
1699     aodPHcells.Sort();
1700   }
1701 }
1702
1703 //______________________________________________________________________________
1704 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)
1705 {
1706   // tracklets
1707   AliCodeTimerAuto("",0);
1708
1709   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());
1710   const AliMultiplicity *mult = esd.GetMultiplicity();
1711   if (mult) {
1712     if (mult->GetNumberOfTracklets()>0) {
1713       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
1714
1715       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
1716         if(fMChandler){
1717           fMChandler->SelectParticle(mult->GetLabel(n, 0));
1718           fMChandler->SelectParticle(mult->GetLabel(n, 1));
1719         }
1720         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1721       }
1722     }
1723   } else {
1724     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
1725   }
1726 }
1727
1728 //______________________________________________________________________________
1729 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)
1730 {
1731   AliCodeTimerAuto("",0);
1732
1733   // Kinks: it is a big mess the access to the information in the kinks
1734   // The loop is on the tracks in order to find the mother and daugther of each kink
1735
1736   Double_t covTr[21]={0.};
1737   Double_t pid[10]={0.};
1738   AliAODPid* detpid(0x0);
1739
1740   fNumberOfKinks = esd.GetNumberOfKinks();
1741
1742   const AliESDVertex* vtx = esd.GetPrimaryVertex();
1743
1744   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
1745   {
1746     AliESDtrack * esdTrack = esd.GetTrack(iTrack);
1747
1748     Int_t ikink = esdTrack->GetKinkIndex(0);
1749
1750     if (ikink && fNumberOfKinks) {
1751             // Negative kink index: mother, positive: daughter
1752
1753             // Search for the second track of the kink
1754
1755             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {
1756
1757         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);
1758
1759         Int_t jkink = esdTrack1->GetKinkIndex(0);
1760
1761         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
1762
1763           // The two tracks are from the same kink
1764
1765           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
1766
1767           Int_t imother = -1;
1768           Int_t idaughter = -1;
1769
1770           if (ikink<0 && jkink>0) {
1771
1772             imother = iTrack;
1773             idaughter = jTrack;
1774           }
1775           else if (ikink>0 && jkink<0) {
1776
1777             imother = jTrack;
1778             idaughter = iTrack;
1779           }
1780           else {
1781             //                  cerr << "Error: Wrong combination of kink indexes: "
1782             //                       << ikink << " " << jkink << endl;
1783             continue;
1784           }
1785
1786           // Add the mother track if it passed primary track selection cuts
1787
1788           AliAODTrack * mother = NULL;
1789
1790           UInt_t selectInfo = 0;
1791           if (fTrackFilter) {
1792             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));
1793             if (!selectInfo) continue;
1794           }
1795
1796           if (!fUsedTrack[imother]) {
1797
1798             fUsedTrack[imother] = kTRUE;
1799
1800             AliESDtrack *esdTrackM = esd.GetTrack(imother);
1801             Double_t p[3] = { 0. };
1802             Double_t pos[3] = { 0. };
1803             esdTrackM->GetPxPyPz(p);
1804             esdTrackM->GetXYZ(pos);
1805             esdTrackM->GetCovarianceXYZPxPyPz(covTr);
1806             esdTrackM->GetESDpid(pid);
1807             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());
1808             mother =
1809             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),
1810                                                esdTrackM->GetLabel(),
1811                                                p,
1812                                                kTRUE,
1813                                                pos,
1814                                                kFALSE,
1815                                                covTr,
1816                                                (Short_t)esdTrackM->GetSign(),
1817                                                esdTrackM->GetITSClusterMap(),
1818                                                pid,
1819                                                fPrimaryVertex,
1820                                                kTRUE, // check if this is right
1821                                                vtx->UsesTrack(esdTrack->GetID()),
1822                                                AliAODTrack::kPrimary,
1823                                                selectInfo);
1824             mother->SetTPCFitMap(esdTrackM->GetTPCFitMap());
1825             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());
1826             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());
1827             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));
1828             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());
1829             mother->SetTPCNCrossedRows(UShort_t(esdTrackM->GetTPCCrossedRows()));
1830
1831             fAODTrackRefs->AddAt(mother, imother);
1832
1833             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;
1834             mother->SetFlags(esdTrackM->GetStatus());
1835             mother->ConvertAliPIDtoAODPID();
1836             fPrimaryVertex->AddDaughter(mother);
1837             mother->ConvertAliPIDtoAODPID();
1838             SetAODPID(esdTrackM,mother,detpid);
1839           }
1840           else {
1841             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1842             //                       << " track " << imother << " has already been used!" << endl;
1843           }
1844
1845           // Add the kink vertex
1846           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);
1847
1848           AliAODVertex * vkink =
1849           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),
1850                                                   NULL,
1851                                                   0.,
1852                                                   mother,
1853                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!
1854                                                   AliAODVertex::kKink);
1855           // Add the daughter track
1856
1857           AliAODTrack * daughter = NULL;
1858
1859           if (!fUsedTrack[idaughter]) {
1860
1861             fUsedTrack[idaughter] = kTRUE;
1862
1863             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);
1864             Double_t p[3] = { 0. };
1865             Double_t pos[3] = { 0. };
1866
1867             esdTrackD->GetPxPyPz(p);
1868             esdTrackD->GetXYZ(pos);
1869             esdTrackD->GetCovarianceXYZPxPyPz(covTr);
1870             esdTrackD->GetESDpid(pid);
1871             selectInfo = 0;
1872             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);
1873             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());
1874             daughter =
1875             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),
1876                                                esdTrackD->GetLabel(),
1877                                                p,
1878                                                kTRUE,
1879                                                pos,
1880                                                kFALSE,
1881                                                covTr,
1882                                                (Short_t)esdTrackD->GetSign(),
1883                                                esdTrackD->GetITSClusterMap(),
1884                                                pid,
1885                                                vkink,
1886                                                kTRUE, // check if this is right
1887                                                vtx->UsesTrack(esdTrack->GetID()),
1888                                                AliAODTrack::kSecondary,
1889                                                selectInfo);
1890             daughter->SetTPCFitMap(esdTrackD->GetTPCFitMap());
1891             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());
1892             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());
1893             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());
1894             daughter->SetTPCNCrossedRows(UShort_t(esdTrackD->GetTPCCrossedRows()));
1895             fAODTrackRefs->AddAt(daughter, idaughter);
1896
1897             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;
1898             daughter->SetFlags(esdTrackD->GetStatus());
1899             daughter->ConvertAliPIDtoAODPID();
1900             vkink->AddDaughter(daughter);
1901             daughter->ConvertAliPIDtoAODPID();
1902             SetAODPID(esdTrackD,daughter,detpid);
1903           }
1904           else {
1905             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1906             //                       << " track " << idaughter << " has already been used!" << endl;
1907           }
1908         }
1909             }
1910     }
1911   }
1912 }
1913
1914 //______________________________________________________________________________
1915 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)
1916 {
1917   AliCodeTimerAuto("",0);
1918
1919   // Access to the AOD container of vertices
1920   fNumberOfVertices = 0;
1921
1922   Double_t pos[3] = { 0. };
1923   Double_t covVtx[6] = { 0. };
1924
1925   // Add primary vertex. The primary tracks will be defined
1926   // after the loops on the composite objects (V0, cascades, kinks)
1927   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1928
1929   vtx->GetXYZ(pos); // position
1930   vtx->GetCovMatrix(covVtx); //covariance matrix
1931
1932   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])
1933   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1934   fPrimaryVertex->SetName(vtx->GetName());
1935   fPrimaryVertex->SetTitle(vtx->GetTitle());
1936   fPrimaryVertex->SetBC(vtx->GetBC());
1937
1938   TString vtitle = vtx->GetTitle();
1939   if (!vtitle.Contains("VertexerTracks"))
1940     fPrimaryVertex->SetNContributors(vtx->GetNContributors());
1941
1942   if (fDebug > 0) fPrimaryVertex->Print();
1943
1944   // Add SPD "main" vertex
1945   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
1946   vtxS->GetXYZ(pos); // position
1947   vtxS->GetCovMatrix(covVtx); //covariance matrix
1948   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])
1949   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
1950   mVSPD->SetName(vtxS->GetName());
1951   mVSPD->SetTitle(vtxS->GetTitle());
1952   mVSPD->SetNContributors(vtxS->GetNContributors());
1953
1954   // Add SPD pileup vertices
1955   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
1956   {
1957     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
1958     vtxP->GetXYZ(pos); // position
1959     vtxP->GetCovMatrix(covVtx); //covariance matrix
1960     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])
1961     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
1962     pVSPD->SetName(vtxP->GetName());
1963     pVSPD->SetTitle(vtxP->GetTitle());
1964     pVSPD->SetNContributors(vtxP->GetNContributors());
1965     pVSPD->SetBC(vtxP->GetBC());
1966   }
1967
1968   // Add TRK pileup vertices
1969   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
1970   {
1971     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
1972     vtxP->GetXYZ(pos); // position
1973     vtxP->GetCovMatrix(covVtx); //covariance matrix
1974     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])
1975     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
1976     pVTRK->SetName(vtxP->GetName());
1977     pVTRK->SetTitle(vtxP->GetTitle());
1978     pVTRK->SetNContributors(vtxP->GetNContributors());
1979     pVTRK->SetBC(vtxP->GetBC());
1980   }
1981
1982   // Add TPC "main" vertex
1983   const AliESDVertex *vtxT = esd.GetPrimaryVertexTPC();
1984   vtxT->GetXYZ(pos); // position
1985   vtxT->GetCovMatrix(covVtx); //covariance matrix
1986   AliAODVertex * mVTPC = new(Vertices()[fNumberOfVertices++])
1987   AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
1988   mVTPC->SetName(vtxT->GetName());
1989   mVTPC->SetTitle(vtxT->GetTitle());
1990   mVTPC->SetNContributors(vtxT->GetNContributors());
1991
1992
1993 }
1994
1995 //______________________________________________________________________________
1996 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)
1997 {
1998   // Convert VZERO data
1999   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
2000   *vzeroData = *(esd.GetVZEROData());
2001 }
2002
2003 //______________________________________________________________________________
2004 void AliAnalysisTaskESDfilter::ConvertTZERO(const AliESDEvent& esd)
2005 {
2006   // Convert TZERO data
2007   const AliESDTZERO* esdTzero = esd.GetESDTZERO();
2008   AliAODTZERO* aodTzero = AODEvent()->GetTZEROData();
2009
2010   for (Int_t icase=0; icase<3; icase++){
2011     aodTzero->SetT0TOF(    icase, esdTzero->GetT0TOF(icase));
2012     aodTzero->SetT0TOFbest(icase, esdTzero->GetT0TOFbest(icase));
2013   }
2014   aodTzero->SetBackgroundFlag(esdTzero->GetBackgroundFlag());
2015   aodTzero->SetPileupFlag(esdTzero->GetPileupFlag());
2016   aodTzero->SetSatelliteFlag(esdTzero->GetSatellite());
2017
2018   Float_t rawTime[24];
2019   for(Int_t ipmt=0; ipmt<24; ipmt++)
2020     rawTime[ipmt] = esdTzero->GetTimeFull(ipmt,0);
2021
2022   Int_t idxOfFirstPmtA = -1,       idxOfFirstPmtC = -1;
2023   Float_t timeOfFirstPmtA = 9999, timeOfFirstPmtC = 9999;
2024   for(int ipmt=0;  ipmt<12; ipmt++){
2025     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtC && rawTime[ipmt]!=0){
2026       timeOfFirstPmtC = rawTime[ipmt];
2027       idxOfFirstPmtC  = ipmt;
2028     }
2029   }
2030   for(int ipmt=12; ipmt<24; ipmt++){
2031     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtA && rawTime[ipmt]!=0 ){
2032       timeOfFirstPmtA = rawTime[ipmt];
2033       idxOfFirstPmtA  = ipmt;
2034     }
2035   }
2036
2037   if(idxOfFirstPmtA != -1 && idxOfFirstPmtC != -1){
2038     //speed of light in cm/ns   TMath::C()*1e-7
2039     Float_t vertexraw = TMath::C()*1e-7 * (rawTime[idxOfFirstPmtA] - rawTime[idxOfFirstPmtC])/2;
2040     aodTzero->SetT0VertexRaw( vertexraw );
2041   }else{
2042     aodTzero->SetT0VertexRaw(99999);
2043   }
2044
2045   aodTzero->SetT0zVertex(esdTzero->GetT0zVertex());
2046 }
2047
2048
2049 //______________________________________________________________________________
2050 void AliAnalysisTaskESDfilter::ConvertZDC(const AliESDEvent& esd)
2051 {
2052   // Convert ZDC data
2053   AliESDZDC* esdZDC = esd.GetZDCData();
2054
2055   const Double_t zem1Energy = esdZDC->GetZEM1Energy();
2056   const Double_t zem2Energy = esdZDC->GetZEM2Energy();
2057
2058   const Double_t *towZNC = esdZDC->GetZNCTowerEnergy();
2059   const Double_t *towZPC = esdZDC->GetZPCTowerEnergy();
2060   const Double_t *towZNA = esdZDC->GetZNATowerEnergy();
2061   const Double_t *towZPA = esdZDC->GetZPATowerEnergy();
2062   const Double_t *towZNCLG = esdZDC->GetZNCTowerEnergyLR();
2063   const Double_t *towZNALG = esdZDC->GetZNATowerEnergyLR();
2064
2065   AliAODZDC* zdcAOD = AODEvent()->GetZDCData();
2066
2067   zdcAOD->SetZEM1Energy(zem1Energy);
2068   zdcAOD->SetZEM2Energy(zem2Energy);
2069   zdcAOD->SetZNCTowers(towZNC, towZNCLG);
2070   zdcAOD->SetZNATowers(towZNA, towZNALG);
2071   zdcAOD->SetZPCTowers(towZPC);
2072   zdcAOD->SetZPATowers(towZPA);
2073
2074   zdcAOD->SetZDCParticipants(esdZDC->GetZDCParticipants(), esdZDC->GetZDCPartSideA(), esdZDC->GetZDCPartSideC());
2075   zdcAOD->SetZDCImpactParameter(esdZDC->GetImpactParameter(), esdZDC->GetImpactParamSideA(),
2076         esdZDC->GetImpactParamSideC());
2077   zdcAOD->SetZDCTDCSum(esdZDC->GetZNTDCSum(0));
2078   zdcAOD->SetZDCTDCDiff(esdZDC->GetZNTDCDiff(0));
2079   if(esdZDC->IsZNChit()) zdcAOD->SetZNCTDC(esdZDC->GetZDCTDCCorrected(10,0));
2080   if(esdZDC->IsZNAhit()) zdcAOD->SetZNATDC(esdZDC->GetZDCTDCCorrected(12,0));
2081 }
2082
2083 //_______________________________________________________________________________________________________________________________________
2084 Int_t AliAnalysisTaskESDfilter::ConvertHMPID(const AliESDEvent& esd) // clm
2085 {
2086   //
2087   // Convtert ESD HMPID info to AOD and return the number of good tracks with HMPID signal.
2088   // We need to return an int since there is no signal counter in the ESD.
2089   //
2090
2091   AliCodeTimerAuto("",0);
2092
2093   Int_t cntHmpidGoodTracks = 0;
2094
2095   Float_t  xMip = 0;
2096   Float_t  yMip = 0;
2097   Int_t    qMip = 0;
2098   Int_t    nphMip = 0;
2099
2100   Float_t xTrk = 0;
2101   Float_t yTrk = 0;
2102   Float_t thetaTrk = 0;
2103   Float_t phiTrk = 0;
2104
2105   Double_t hmpPid[5]={0};
2106   Double_t hmpMom[3]={0};
2107
2108   TClonesArray &hmpidRings = *(AODEvent()->GetHMPIDrings());
2109
2110  for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
2111   {
2112     if(! esd.GetTrack(iTrack) ) continue;
2113
2114     if(esd.GetTrack(iTrack)->GetHMPIDsignal() > -20 ) {                  //
2115
2116       (esd.GetTrack(iTrack))->GetHMPIDmip(xMip, yMip, qMip, nphMip);     // Get MIP properties
2117       (esd.GetTrack(iTrack))->GetHMPIDtrk(xTrk,yTrk,thetaTrk,phiTrk);
2118       (esd.GetTrack(iTrack))->GetHMPIDpid(hmpPid);
2119       if((esd.GetTrack(iTrack))->GetOuterHmpParam()) (esd.GetTrack(iTrack))->GetOuterHmpPxPyPz(hmpMom);
2120
2121      if(esd.GetTrack(iTrack)->GetHMPIDsignal() == 0 && thetaTrk == 0 && qMip == 0 && nphMip ==0 ) continue;  //
2122
2123      new(hmpidRings[cntHmpidGoodTracks++]) AliAODHMPIDrings(
2124                                                                  (esd.GetTrack(iTrack))->GetID(),             // Unique track id to attach the ring to
2125                                                                  1000000*nphMip+qMip,                         // MIP charge and number of photons
2126                                                                  (esd.GetTrack(iTrack))->GetHMPIDcluIdx(),    // 1000000*chamber id + cluster idx of the assigned MIP cluster
2127                                                                  thetaTrk,                                    // track inclination angle theta
2128                                                                  phiTrk,                                      // track inclination angle phi
2129                                                                  (esd.GetTrack(iTrack))->GetHMPIDsignal(),    // Cherenkov angle
2130                                                                  (esd.GetTrack(iTrack))->GetHMPIDoccupancy(), // Occupancy claculated for the given chamber
2131                                                                  (esd.GetTrack(iTrack))->GetHMPIDchi2(),      // Ring resolution squared
2132                                                                   xTrk,                                       // Track x coordinate (LORS)
2133                                                                   yTrk,                                       // Track y coordinate (LORS)
2134                                                                   xMip,                                       // MIP x coordinate (LORS)
2135                                                                   yMip,                                       // MIP y coordinate (LORS)
2136                                                                   hmpPid,                                     // PID probablities from ESD, remove later once it is in CombinedPid
2137                                                                   hmpMom                                      // Track momentum in HMPID at ring reconstruction
2138                                                                );
2139
2140       //  Printf(Form("+++++++++ yes/no: %d  %lf %lf %lf %lf %lf %lf ",(esd.GetTrack(iTrack))->IsHMPID(),thetaTrk, (esd.GetTrack(iTrack))->GetHMPIDchi2(),xTrk, yTrk , xMip, yMip));
2141
2142
2143    }// HMPID signal > -20
2144   }//___esd track loop
2145
2146   return cntHmpidGoodTracks;
2147 }
2148
2149 //______________________________________________________________________________
2150 void AliAnalysisTaskESDfilter::ConvertESDtoAOD()
2151 {
2152   // ESD Filter analysis task executed for each event
2153
2154   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
2155
2156   if(!esd)return;
2157
2158   AliCodeTimerAuto("",0);
2159
2160   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );
2161
2162       // Reconstruct cascades and V0 here
2163   if (fIsV0CascadeRecoEnabled) {
2164     esd->ResetCascades();
2165     esd->ResetV0s();
2166
2167     AliV0vertexer lV0vtxer;
2168     AliCascadeVertexer lCascVtxer;
2169
2170     lV0vtxer.SetCuts(fV0Cuts);
2171     lCascVtxer.SetCuts(fCascadeCuts);
2172
2173
2174     lV0vtxer.Tracks2V0vertices(esd);
2175     lCascVtxer.V0sTracks2CascadeVertices(esd);
2176   }
2177
2178
2179   fNumberOfTracks = 0;
2180   fNumberOfPositiveTracks = 0;
2181   fNumberOfV0s = 0;
2182   fNumberOfVertices = 0;
2183   fNumberOfCascades = 0;
2184   fNumberOfKinks = 0;
2185
2186   AliAODHeader* header = ConvertHeader(*esd);
2187
2188   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);
2189   if ( fIsTZEROEnabled ) ConvertTZERO(*esd);
2190
2191   // Fetch Stack for debuggging if available
2192   fMChandler=0x0;
2193   if(MCEvent())
2194   {
2195     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
2196   }
2197
2198   // loop over events and fill them
2199   // Multiplicity information needed by the header (to be revised!)
2200   Int_t nTracks    = esd->GetNumberOfTracks();
2201   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);
2202
2203   // Update the header
2204
2205   Int_t nV0s      = esd->GetNumberOfV0s();
2206   Int_t nCascades = esd->GetNumberOfCascades();
2207   Int_t nKinks    = esd->GetNumberOfKinks();
2208   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
2209   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
2210   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
2211   nVertices+=nPileSPDVertices;
2212   nVertices+=nPileTrkVertices;
2213   Int_t nJets     = 0;
2214   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
2215   Int_t nFmdClus  = 0;
2216   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
2217   Int_t nHmpidRings = 0;
2218
2219   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));
2220
2221   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus,nHmpidRings);
2222
2223   if (nV0s > 0)
2224   {
2225     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0
2226     fAODV0VtxRefs = new TRefArray(nV0s);
2227     // RefArray to store the mapping between esd V0 number and newly created AOD-V0
2228     fAODV0Refs = new TRefArray(nV0s);
2229     // Array to take into account the V0s already added to the AOD (V0 within cascades)
2230     fUsedV0 = new Bool_t[nV0s];
2231     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;
2232   }
2233
2234   if (nTracks>0)
2235   {
2236     // RefArray to store the mapping between esd track number and newly created AOD-Track
2237
2238     fAODTrackRefs = new TRefArray(nTracks);
2239
2240     // Array to take into account the tracks already added to the AOD
2241     fUsedTrack = new Bool_t[nTracks];
2242     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;
2243   }
2244
2245   // Array to take into account the kinks already added to the AOD
2246   if (nKinks>0)
2247   {
2248     fUsedKink = new Bool_t[nKinks];
2249     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;
2250   }
2251
2252   ConvertPrimaryVertices(*esd);
2253
2254   //setting best TOF PID
2255   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
2256   if (esdH)
2257       fESDpid = esdH->GetESDpid();
2258
2259   if (fIsPidOwner && fESDpid){
2260     delete fESDpid;
2261     fESDpid = 0;
2262   }
2263   if(!fESDpid)
2264   { //in case of no Tender attached
2265     fESDpid = new AliESDpid;
2266     fIsPidOwner = kTRUE;
2267   }
2268
2269   if(!esd->GetTOFHeader())
2270   { //protection in case the pass2 LHC10b,c,d have been processed without tender.
2271     Float_t t0spread[10];
2272     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!!
2273     for (Int_t i=0; i<10; i++) t0spread[i] = (TMath::Sqrt(esd->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
2274     fESDpid->GetTOFResponse().SetT0resolution(t0spread);
2275     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
2276       //    fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);
2277     AliTOFHeader tmpTOFHeader(0,t0spread[0],0,NULL,NULL,NULL,intrinsicTOFres,t0spread[0]);
2278     AODEvent()->SetTOFHeader(&tmpTOFHeader);         // write dummy TOF header in AOD
2279   } else {
2280     AODEvent()->SetTOFHeader(esd->GetTOFHeader());    // write TOF header in AOD
2281   }
2282
2283   //  if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender.
2284
2285   if ( fAreCascadesEnabled ) ConvertCascades(*esd);
2286
2287   if ( fAreV0sEnabled ) ConvertV0s(*esd);
2288
2289   if ( fAreKinksEnabled ) ConvertKinks(*esd);
2290
2291   if ( fAreTracksEnabled ) ConvertTracks(*esd);
2292
2293   // Update number of AOD tracks in header at the end of track loop (M.G.)
2294   header->SetRefMultiplicity(fNumberOfTracks);
2295   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);
2296   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);
2297
2298   if ( fTPCConstrainedFilterMask ) ConvertTPCOnlyTracks(*esd);
2299   if( fGlobalConstrainedFilterMask) ConvertGlobalConstrainedTracks(*esd);
2300
2301   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);
2302
2303   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);
2304
2305   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);
2306
2307   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);
2308
2309         if ( fAreEMCALTriggerEnabled )ConvertCaloTrigger(TString("EMCAL"), *esd);
2310
2311         if ( fArePHOSTriggerEnabled )ConvertCaloTrigger(TString("PHOS"), *esd);
2312
2313   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);
2314   if ( fIsZDCEnabled ) ConvertZDC(*esd);
2315
2316   if(fIsHMPIDEnabled) nHmpidRings = ConvertHMPID(*esd);
2317
2318   delete fAODTrackRefs; fAODTrackRefs=0x0;
2319   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;
2320   delete fAODV0Refs; fAODV0Refs=0x0;
2321
2322   delete[] fUsedTrack; fUsedTrack=0x0;
2323   delete[] fUsedV0; fUsedV0=0x0;
2324   delete[] fUsedKink; fUsedKink=0x0;
2325
2326   if ( fIsPidOwner){
2327     delete fESDpid;
2328     fESDpid = 0x0;
2329   }
2330
2331
2332 }
2333
2334
2335 //______________________________________________________________________________
2336 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)
2337 {
2338   //
2339   // Setter for the raw PID detector signals
2340   //
2341
2342   // Save PID object for candidate electrons
2343     Bool_t pidSave = kFALSE;
2344     if (fTrackFilter) {
2345         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");
2346         if (selectInfo)  pidSave = kTRUE;
2347     }
2348
2349
2350     // Tracks passing pt cut
2351     if(esdtrack->Pt()>fHighPthreshold) {
2352         pidSave = kTRUE;
2353     } else {
2354         if(fPtshape){
2355             if(esdtrack->Pt()> fPtshape->GetXmin()){
2356                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
2357                 if(gRandom->Rndm(0)<1./y){
2358                     pidSave = kTRUE;
2359                 }//end rndm
2360             }//end if p < pmin
2361         }//end if p function
2362     }// end else
2363
2364     if (pidSave) {
2365       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track
2366         detpid = new AliAODPid();
2367         SetDetectorRawSignals(detpid,esdtrack);
2368         aodtrack->SetDetPID(detpid);
2369       }
2370     }
2371 }
2372
2373 //______________________________________________________________________________
2374 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)
2375 {
2376 //
2377 //assignment of the detector signals (AliXXXesdPID inspired)
2378 //
2379  if(!track){
2380  AliInfo("no ESD track found. .....exiting");
2381  return;
2382  }
2383  // TPC momentum
2384  const AliExternalTrackParam *in=track->GetInnerParam();
2385  if (in) {
2386    aodpid->SetTPCmomentum(in->GetP());
2387  }else{
2388    aodpid->SetTPCmomentum(-1.);
2389  }
2390
2391
2392  aodpid->SetITSsignal(track->GetITSsignal());
2393  Double_t itsdedx[4]; // dE/dx samples for individual ITS layers
2394  track->GetITSdEdxSamples(itsdedx);
2395  aodpid->SetITSdEdxSamples(itsdedx);
2396
2397  aodpid->SetTPCsignal(track->GetTPCsignal());
2398  aodpid->SetTPCsignalN(track->GetTPCsignalN());
2399  if(track->GetTPCdEdxInfo()) aodpid->SetTPCdEdxInfo(track->GetTPCdEdxInfo());
2400
2401  //n TRD planes = 6
2402  Int_t nslices = track->GetNumberOfTRDslices()*6;
2403  TArrayD trdslices(nslices);
2404  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
2405    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
2406  }
2407
2408 //TRD momentum
2409  for(Int_t iPl=0;iPl<6;iPl++){
2410    Double_t trdmom=track->GetTRDmomentum(iPl);
2411    aodpid->SetTRDmomentum(iPl,trdmom);
2412  }
2413
2414  aodpid->SetTRDslices(track->GetNumberOfTRDslices()*6,trdslices.GetArray());
2415  aodpid->SetTRDsignal(track->GetTRDsignal());
2416
2417  //TRD clusters and tracklets
2418  aodpid->SetTRDncls(track->GetTRDncls());
2419  aodpid->SetTRDntrackletsPID(track->GetTRDntrackletsPID());
2420
2421  aodpid->SetTRDChi2(track->GetTRDchi2());
2422
2423  //TOF PID
2424  Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
2425  aodpid->SetIntegratedTimes(times);
2426
2427    //  Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
2428    //  aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);
2429    aodpid->SetTOFsignal(track->GetTOFsignal());
2430
2431   Double_t tofRes[5];
2432   for (Int_t iMass=0; iMass<5; iMass++){
2433     //    tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));
2434     tofRes[iMass]=0; //backward compatibility
2435   }
2436   aodpid->SetTOFpidResolution(tofRes);
2437
2438 //  aodpid->SetHMPIDsignal(0); // set to zero for compression but it will be removed later
2439
2440 }
2441
2442 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)
2443 {
2444     // Calculate chi2 per ndf for track
2445     Int_t  nClustersTPC = track->GetTPCNcls();
2446
2447     if ( nClustersTPC > 5) {
2448        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));
2449     } else {
2450        return (-1.);
2451     }
2452  }
2453
2454
2455 //______________________________________________________________________________
2456 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
2457 {
2458 // Terminate analysis
2459 //
2460     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
2461 }
2462
2463 //______________________________________________________________________________
2464 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){
2465 // Print MC info
2466   if(!pStack)return;
2467   label = TMath::Abs(label);
2468   TParticle *part = pStack->Particle(label);
2469   Printf("########################");
2470   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
2471   part->Print();
2472   TParticle* mother = part;
2473   Int_t imo = part->GetFirstMother();
2474   Int_t nprim = pStack->GetNprimary();
2475   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
2476   while((imo >= nprim)) {
2477     mother =  pStack->Particle(imo);
2478     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
2479     mother->Print();
2480     imo =  mother->GetFirstMother();
2481   }
2482   Printf("########################");
2483 }
2484
2485 //______________________________________________________
2486