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