ANALYSIS/AliAnalysisTaskESDfilter.cxx

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