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