ANALYSIS/AliAnalysisTaskESDfilter.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 /* $Id$ */\r
  17  \r
  18 #include <TChain.h>\r
  19 #include <TTree.h>\r
  20 #include <TList.h>\r
  21 #include <TArrayI.h>\r
  22 #include <TRandom.h>\r
  23 #include <TParticle.h>\r
  24 #include <TFile.h>\r
  25 \r
  26 #include "AliAnalysisTaskESDfilter.h"\r
  27 #include "AliAnalysisManager.h"\r
  28 #include "AliESDEvent.h"\r
  29 #include "AliESDRun.h"\r
  30 #include "AliStack.h"\r
  31 #include "AliAODEvent.h"\r
  32 #include "AliMCEvent.h"\r
  33 #include "AliMCEventHandler.h"\r
  34 #include "AliESDInputHandler.h"\r
  35 #include "AliAODHandler.h"\r
  36 #include "AliAODMCParticle.h"\r
  37 #include "AliAnalysisFilter.h"\r
  38 #include "AliESDMuonTrack.h"\r
  39 #include "AliESDVertex.h"\r
  40 #include "AliCentrality.h"\r
  41 #include "AliEventplane.h"\r
  42 #include "AliESDv0.h"\r
  43 #include "AliESDkink.h"\r
  44 #include "AliESDcascade.h"\r
  45 #include "AliESDPmdTrack.h"\r
  46 #include "AliESDCaloCluster.h"\r
  47 #include "AliESDCaloCells.h"\r
  48 #include "AliMultiplicity.h"\r
  49 #include "AliLog.h"\r
  50 #include "AliCodeTimer.h"\r
  51 #include "AliESDtrackCuts.h"\r
  52 #include "AliESDpid.h"\r
  53 #include "Riostream.h"\r
  54 \r
  55 ClassImp(AliAnalysisTaskESDfilter)\r
  56 \r
  57 ////////////////////////////////////////////////////////////////////////\r
  58 \r
  59 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():\r
  60     AliAnalysisTaskSE(),\r
  61     fTrackFilter(0x0),\r
  62     fKinkFilter(0x0),\r
  63     fV0Filter(0x0),\r
  64     fCascadeFilter(0x0),\r
  65     fHighPthreshold(0),\r
  66     fPtshape(0x0),\r
  67     fEnableFillAOD(kTRUE),\r
  68 fUsedTrack(0x0),\r
  69 fUsedKink(0x0),\r
  70 fUsedV0(0x0),\r
  71 fAODTrackRefs(0x0),\r
  72 fAODV0VtxRefs(0x0),\r
  73 fAODV0Refs(0x0),\r
  74 fMChandler(0x0),\r
  75 fNumberOfTracks(0),\r
  76 fNumberOfPositiveTracks(0),\r
  77 fNumberOfV0s(0),\r
  78 fNumberOfVertices(0),\r
  79 fNumberOfCascades(0),\r
  80 fNumberOfKinks(0),\r
  81 fOldESDformat(kFALSE),\r
  82 fPrimaryVertex(0x0),\r
  83 fTPCOnlyFilterMask(0),\r
  84 fIsVZEROEnabled(kTRUE),\r
  85 fAreCascadesEnabled(kTRUE),\r
  86 fAreV0sEnabled(kTRUE),\r
  87 fAreKinksEnabled(kTRUE),\r
  88 fAreTracksEnabled(kTRUE),\r
  89 fArePmdClustersEnabled(kTRUE),\r
  90 fAreCaloClustersEnabled(kTRUE),\r
  91 fAreEMCALCellsEnabled(kTRUE),\r
  92 fArePHOSCellsEnabled(kTRUE),\r
  93 fAreTrackletsEnabled(kTRUE),\r
  94 fESDpid(0x0),\r
  95 fIsPidOwner(kFALSE),\r
  96 fTimeZeroType(AliESDpid::kTOF_T0)\r
  97 {\r
  98   // Default constructor\r
  99 }\r
 100 \r
 101 //______________________________________________________________________________\r
 102 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):\r
 103     AliAnalysisTaskSE(name),\r
 104     fTrackFilter(0x0),\r
 105     fKinkFilter(0x0),\r
 106     fV0Filter(0x0),\r
 107     fCascadeFilter(0x0),\r
 108     fHighPthreshold(0),\r
 109     fPtshape(0x0),\r
 110     fEnableFillAOD(kTRUE),\r
 111 fUsedTrack(0x0),\r
 112 fUsedKink(0x0),\r
 113 fUsedV0(0x0),\r
 114 fAODTrackRefs(0x0),\r
 115 fAODV0VtxRefs(0x0),\r
 116 fAODV0Refs(0x0),\r
 117 fMChandler(0x0),\r
 118 fNumberOfTracks(0),\r
 119 fNumberOfPositiveTracks(0),\r
 120 fNumberOfV0s(0),\r
 121 fNumberOfVertices(0),\r
 122 fNumberOfCascades(0),\r
 123 fNumberOfKinks(0),\r
 124 fOldESDformat(kFALSE),\r
 125 fPrimaryVertex(0x0),\r
 126 fTPCOnlyFilterMask(0),\r
 127 fIsVZEROEnabled(kTRUE),\r
 128 fAreCascadesEnabled(kTRUE),\r
 129 fAreV0sEnabled(kTRUE),\r
 130 fAreKinksEnabled(kTRUE),\r
 131 fAreTracksEnabled(kTRUE),\r
 132 fArePmdClustersEnabled(kTRUE),\r
 133 fAreCaloClustersEnabled(kTRUE),\r
 134 fAreEMCALCellsEnabled(kTRUE),\r
 135 fArePHOSCellsEnabled(kTRUE),\r
 136 fAreTrackletsEnabled(kTRUE),\r
 137 fESDpid(0x0),\r
 138 fIsPidOwner(kFALSE),\r
 139 fTimeZeroType(AliESDpid::kTOF_T0)\r
 140 {\r
 141   // Constructor\r
 142 }\r
 143 AliAnalysisTaskESDfilter::~AliAnalysisTaskESDfilter(){\r
 144   if(fIsPidOwner)delete fESDpid;\r
 145 }\r
 146 //______________________________________________________________________________\r
 147 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()\r
 148 {\r
 149   //\r
 150   // Create Output Objects conenct filter to outputtree\r
 151   // \r
 152   if(OutputTree())\r
 153   {\r
 154     OutputTree()->GetUserInfo()->Add(fTrackFilter);\r
 155   }\r
 156   else\r
 157   {\r
 158     AliError("No OutputTree() for adding the track filter");\r
 159   }\r
 160 }\r
 161 \r
 162 //______________________________________________________________________________\r
 163 void AliAnalysisTaskESDfilter::Init()\r
 164 {\r
 165   // Initialization\r
 166   if (fDebug > 1) AliInfo("Init() \n");\r
 167   // Call configuration file\r
 168 }\r
 169 \r
 170 //______________________________________________________________________________\r
 171 void AliAnalysisTaskESDfilter::PrintTask(Option_t *option, Int_t indent) const\r
 172 {\r
 173   AliInfo("");\r
 174   \r
 175   AliAnalysisTaskSE::PrintTask(option,indent);\r
 176   \r
 177   TString spaces(' ',indent+3);\r
 178   \r
 179   cout << spaces.Data() << Form("Cascades     are %s",fAreCascadesEnabled ? "ENABLED":"DISABLED") << endl;\r
 180   cout << spaces.Data() << Form("V0s          are %s",fAreV0sEnabled ? "ENABLED":"DISABLED") << endl;\r
 181   cout << spaces.Data() << Form("Kinks        are %s",fAreKinksEnabled ? "ENABLED":"DISABLED") << endl;\r
 182   cout << spaces.Data() << Form("Tracks       are %s",fAreTracksEnabled ? "ENABLED":"DISABLED") << endl;\r
 183   cout << spaces.Data() << Form("PmdClusters  are %s",fArePmdClustersEnabled ? "ENABLED":"DISABLED") << endl;\r
 184   cout << spaces.Data() << Form("CaloClusters are %s",fAreCaloClustersEnabled ? "ENABLED":"DISABLED") << endl;\r
 185   cout << spaces.Data() << Form("EMCAL cells  are %s",fAreEMCALCellsEnabled ? "ENABLED":"DISABLED") << endl;\r
 186   cout << spaces.Data() << Form("Tracklets    are %s",fAreTrackletsEnabled ? "ENABLED":"DISABLED") << endl;  \r
 187 }\r
 188 \r
 189 //______________________________________________________________________________\r
 190 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)\r
 191 {\r
 192 // Execute analysis for current event\r
 193 //\r
 194                                             \r
 195   Long64_t ientry = Entry();\r
 196   \r
 197   if (fDebug > 0) {\r
 198       printf("Filter: Analysing event # %5d\n", (Int_t) ientry);\r
 199       if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");\r
 200       if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");\r
 201   }\r
 202   // Filters must explicitely enable AOD filling in their UserExec (AG)\r
 203   if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");\r
 204   if(fEnableFillAOD) {\r
 205      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);\r
 206      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);\r
 207   }   \r
 208   ConvertESDtoAOD();\r
 209 }\r
 210 \r
 211 //______________________________________________________________________________\r
 212 TClonesArray& AliAnalysisTaskESDfilter::Cascades()\r
 213 {\r
 214   return *(AODEvent()->GetCascades());\r
 215 }\r
 216 \r
 217 //______________________________________________________________________________\r
 218 TClonesArray& AliAnalysisTaskESDfilter::Tracks()\r
 219 {\r
 220   return *(AODEvent()->GetTracks());\r
 221 }\r
 222 \r
 223 //______________________________________________________________________________\r
 224 TClonesArray& AliAnalysisTaskESDfilter::V0s()\r
 225 {\r
 226   return *(AODEvent()->GetV0s());\r
 227 }\r
 228 \r
 229 //______________________________________________________________________________\r
 230 TClonesArray& AliAnalysisTaskESDfilter::Vertices()\r
 231 {\r
 232   return *(AODEvent()->GetVertices());\r
 233 }\r
 234 \r
 235 //______________________________________________________________________________\r
 236 AliAODHeader* AliAnalysisTaskESDfilter::ConvertHeader(const AliESDEvent& esd)\r
 237 {\r
 238   AliCodeTimerAuto("",0);\r
 239   \r
 240   AliAODHeader* header = AODEvent()->GetHeader();\r
 241   \r
 242   header->SetRunNumber(esd.GetRunNumber());\r
 243   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection\r
 244 \r
 245   TTree* tree = fInputHandler->GetTree();\r
 246   if (tree) {\r
 247     TFile* file = tree->GetCurrentFile();\r
 248     if (file) header->SetESDFileName(file->GetName());\r
 249   }\r
 250   \r
 251   if (fOldESDformat) {\r
 252     header->SetBunchCrossNumber(0);\r
 253     header->SetOrbitNumber(0);\r
 254     header->SetPeriodNumber(0);\r
 255     header->SetEventType(0);\r
 256     header->SetMuonMagFieldScale(-999.);\r
 257     header->SetCentrality(0);       \r
 258     header->SetEventplane(0);\r
 259   } else {\r
 260     header->SetBunchCrossNumber(esd.GetBunchCrossNumber());\r
 261     header->SetOrbitNumber(esd.GetOrbitNumber());\r
 262     header->SetPeriodNumber(esd.GetPeriodNumber());\r
 263     header->SetEventType(esd.GetEventType());\r
 264     \r
 265     header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());\r
 266     if(const_cast<AliESDEvent&>(esd).GetCentrality()){\r
 267       header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());\r
 268     }\r
 269     else{\r
 270       header->SetCentrality(0);\r
 271     }\r
 272     if(const_cast<AliESDEvent&>(esd).GetEventplane()){\r
 273       header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());\r
 274     }\r
 275     else{\r
 276       header->SetEventplane(0);\r
 277     }\r
 278   }\r
 279   \r
 280   // Trigger\r
 281   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());\r
 282   header->SetTriggerMask(esd.GetTriggerMask()); \r
 283   header->SetTriggerCluster(esd.GetTriggerCluster());\r
 284   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());    \r
 285   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());    \r
 286   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());    \r
 287   \r
 288   header->SetMagneticField(esd.GetMagneticField());\r
 289   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);\r
 290   header->SetZDCN1Energy(esd.GetZDCN1Energy());\r
 291   header->SetZDCP1Energy(esd.GetZDCP1Energy());\r
 292   header->SetZDCN2Energy(esd.GetZDCN2Energy());\r
 293   header->SetZDCP2Energy(esd.GetZDCP2Energy());\r
 294   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));\r
 295   \r
 296   // ITS Cluster Multiplicty\r
 297   const AliMultiplicity *mult = esd.GetMultiplicity();\r
 298   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));\r
 299   \r
 300   // TPC only Reference Multiplicty\r
 301   AliESDtrackCuts* cuts  = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();\r
 302   Int_t refMult  = cuts ? (Short_t)cuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;\r
 303   header->SetTPConlyRefMultiplicity(refMult);\r
 304   \r
 305   //\r
 306   Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};\r
 307   Float_t diamcov[3]; \r
 308   esd.GetDiamondCovXY(diamcov);\r
 309   header->SetDiamond(diamxy,diamcov);\r
 310   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());\r
 311   \r
 312   return header;\r
 313 }\r
 314 \r
 315 //______________________________________________________________________________\r
 316 void AliAnalysisTaskESDfilter::ConvertCascades(const AliESDEvent& esd) \r
 317 {\r
 318   // Convert the cascades part of the ESD.\r
 319   // Return the number of cascades\r
 320  \r
 321   AliCodeTimerAuto("",0);\r
 322   \r
 323   // Create vertices starting from the most complex objects\r
 324   Double_t chi2 = 0.;\r
 325   \r
 326   const AliESDVertex* vtx = esd.GetPrimaryVertex();\r
 327   Double_t pos[3] = { 0. };\r
 328   Double_t covVtx[6] = { 0. };\r
 329   Double_t momBach[3]={0.};\r
 330   Double_t covTr[21]={0.};\r
 331   Double_t pid[10]={0.};\r
 332   AliAODPid* detpid(0x0);\r
 333   AliAODVertex* vV0FromCascade(0x0);\r
 334   AliAODv0* aodV0(0x0);\r
 335   AliAODcascade* aodCascade(0x0);\r
 336   AliAODTrack* aodTrack(0x0);\r
 337   Double_t momPos[3]={0.};\r
 338   Double_t momNeg[3] = { 0. };\r
 339   Double_t momPosAtV0vtx[3]={0.};\r
 340   Double_t momNegAtV0vtx[3]={0.};\r
 341 \r
 342   TClonesArray& verticesArray = Vertices();\r
 343   TClonesArray& tracksArray = Tracks();\r
 344   TClonesArray& cascadesArray = Cascades();\r
 345   \r
 346   // Cascades (Modified by A.Maire - February 2009)\r
 347   for (Int_t nCascade = 0; nCascade < esd.GetNumberOfCascades(); ++nCascade) {\r
 348     \r
 349     // 0- Preparation\r
 350     //\r
 351     AliESDcascade *esdCascade = esd.GetCascade(nCascade);\r
 352                 Int_t  idxPosFromV0Dghter  = esdCascade->GetPindex();\r
 353                 Int_t  idxNegFromV0Dghter  = esdCascade->GetNindex();\r
 354                 Int_t  idxBachFromCascade  = esdCascade->GetBindex();\r
 355     \r
 356     AliESDtrack  *esdCascadePos  = esd.GetTrack( idxPosFromV0Dghter);\r
 357     AliESDtrack  *esdCascadeNeg  = esd.GetTrack( idxNegFromV0Dghter);\r
 358     AliESDtrack  *esdCascadeBach = esd.GetTrack( idxBachFromCascade);\r
 359     \r
 360     // Identification of the V0 within the esdCascade (via both daughter track indices)\r
 361     AliESDv0 * currentV0   = 0x0;\r
 362     Int_t      idxV0FromCascade = -1;\r
 363     \r
 364     for (Int_t iV0=0; iV0<esd.GetNumberOfV0s(); ++iV0) {\r
 365       \r
 366       currentV0 = esd.GetV0(iV0);\r
 367       Int_t posCurrentV0 = currentV0->GetPindex();\r
 368       Int_t negCurrentV0 = currentV0->GetNindex();\r
 369       \r
 370       if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {\r
 371         idxV0FromCascade = iV0;\r
 372         break;\r
 373       }\r
 374     }\r
 375     \r
 376     if(idxV0FromCascade < 0){\r
 377       printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");\r
 378       continue;\r
 379     }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"\r
 380     \r
 381     AliESDv0 *esdV0FromCascade   = esd.GetV0(idxV0FromCascade);\r
 382         \r
 383     // 1 - Cascade selection \r
 384     \r
 385     //  AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd.GetPrimaryVertex()));\r
 386     //  TList cascadeObjects;\r
 387     //  cascadeObjects.AddAt(esdV0FromCascade, 0);\r
 388     //  cascadeObjects.AddAt(esdCascadePos,    1);\r
 389     //  cascadeObjects.AddAt(esdCascadeNeg,    2);\r
 390     //  cascadeObjects.AddAt(esdCascade,       3);\r
 391     //  cascadeObjects.AddAt(esdCascadeBach,   4);\r
 392     //  cascadeObjects.AddAt(esdPrimVtx,       5);\r
 393     // \r
 394     //  UInt_t selectCascade = 0;\r
 395     //  if (fCascadeFilter) {\r
 396     //    // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects); \r
 397     //          // FIXME AliESDCascadeCuts to be implemented ...\r
 398     // \r
 399     //          // Here we may encounter a moot point at the V0 level \r
 400     //          // between the cascade selections and the V0 ones :\r
 401     //          // the V0 selected along with the cascade (secondary V0) may \r
 402     //          // usually be removed from the dedicated V0 selections (prim V0) ...\r
 403     //          // -> To be discussed !\r
 404     // \r
 405     //    // this is a little awkward but otherwise the \r
 406     //    // list wants to access the pointer (delete it) \r
 407     //    // again when going out of scope\r
 408     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct\r
 409     //    esdPrimVtx = 0;\r
 410     //    if (!selectCascade) \r
 411     //      continue;\r
 412     //  }\r
 413     //  else{\r
 414     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct\r
 415     //    esdPrimVtx = 0;\r
 416     //  }\r
 417     \r
 418     // 2 - Add the cascade vertex\r
 419     \r
 420     esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);\r
 421     esdCascade->GetPosCovXi(covVtx);\r
 422     chi2 = esdCascade->GetChi2Xi(); \r
 423     \r
 424     AliAODVertex *vCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex( pos,\r
 425                                                                      covVtx,\r
 426                                                                      chi2, // FIXME = Chi2/NDF will be needed\r
 427                                                                      fPrimaryVertex,\r
 428                                                                      nCascade, // id\r
 429                                                                      AliAODVertex::kCascade);\r
 430     fPrimaryVertex->AddDaughter(vCascade);\r
 431     \r
 432 //    if (fDebug > 2) {\r
 433 //      printf("---- Cascade / Cascade Vertex (AOD) : \n");\r
 434 //      vCascade->Print();\r
 435 //    }\r
 436     \r
 437     if(esd.GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(const_cast<AliESDEvent*>(&esd), AliESDpid::kBest_T0); //in case of AOD production strating form LHC10e without Tender. \r
 438 \r
 439 \r
 440     // 3 - Add the bachelor track from the cascade\r
 441     \r
 442     if (!fUsedTrack[idxBachFromCascade]) {\r
 443       \r
 444       esdCascadeBach->GetPxPyPz(momBach);\r
 445       esdCascadeBach->GetXYZ(pos);\r
 446       esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);\r
 447       esdCascadeBach->GetESDpid(pid);\r
 448       \r
 449             fUsedTrack[idxBachFromCascade] = kTRUE;\r
 450             UInt_t selectInfo = 0;\r
 451             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);\r
 452             if (fMChandler) fMChandler->SelectParticle(esdCascadeBach->GetLabel());\r
 453             aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(esdCascadeBach->GetID(),\r
 454                                                                            esdCascadeBach->GetLabel(), \r
 455                                                                            momBach, \r
 456                                                                            kTRUE,\r
 457                                                                            pos,\r
 458                                                                            kFALSE, // Why kFALSE for "isDCA" ? FIXME\r
 459                                                                            covTr, \r
 460                                                                            (Short_t)esdCascadeBach->GetSign(),\r
 461                                                                            esdCascadeBach->GetITSClusterMap(), \r
 462                                                                            pid,\r
 463                                                                            vCascade,\r
 464                                                                            kTRUE,  // usedForVtxFit = kFALSE ? FIXME\r
 465                                                                            vtx->UsesTrack(esdCascadeBach->GetID()),\r
 466                                                                            AliAODTrack::kSecondary,\r
 467                                                                            selectInfo);\r
 468             aodTrack->SetTPCClusterMap(esdCascadeBach->GetTPCClusterMap());\r
 469             aodTrack->SetTPCSharedMap (esdCascadeBach->GetTPCSharedMap());\r
 470             aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeBach));\r
 471             aodTrack->SetTPCPointsF(esdCascadeBach->GetTPCNclsF());\r
 472             fAODTrackRefs->AddAt(aodTrack,idxBachFromCascade);\r
 473             \r
 474             if (esdCascadeBach->GetSign() > 0) ++fNumberOfPositiveTracks;\r
 475             aodTrack->ConvertAliPIDtoAODPID();\r
 476             aodTrack->SetFlags(esdCascadeBach->GetStatus());\r
 477       SetAODPID(esdCascadeBach,aodTrack,detpid);\r
 478     }\r
 479     else {\r
 480             aodTrack = static_cast<AliAODTrack*>( fAODTrackRefs->At(idxBachFromCascade) );\r
 481     }\r
 482     \r
 483     vCascade->AddDaughter(aodTrack);\r
 484     \r
 485 //    if (fDebug > 4) {\r
 486 //      printf("---- Cascade / bach dghter : \n");\r
 487 //      aodTrack->Print();\r
 488 //    }\r
 489     \r
 490     \r
 491     // 4 - Add the V0 from the cascade. \r
 492     // = V0vtx + both pos and neg daughter tracks + the aodV0 itself\r
 493     //\r
 494     \r
 495     if ( !fUsedV0[idxV0FromCascade] ) {\r
 496       // 4.A - if VO structure hasn't been created yet\r
 497       \r
 498       // 4.A.1 - Create the V0 vertex of the cascade\r
 499       \r
 500       esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);\r
 501       esdV0FromCascade->GetPosCov(covVtx);\r
 502       chi2 = esdV0FromCascade->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3 ?\r
 503                         \r
 504       vV0FromCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex(pos,\r
 505                                                                          covVtx,\r
 506                                                                          chi2,\r
 507                                                                          vCascade,\r
 508                                                                          idxV0FromCascade, //id of ESDv0\r
 509                                                                          AliAODVertex::kV0);\r
 510       // Note:\r
 511       //    one V0 can be used by several cascades.\r
 512       // So, one AOD V0 vtx can have several parent vtx.\r
 513       // This is not directly allowed by AliAODvertex.\r
 514       // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash\r
 515       // but to a problem of consistency within AODEvent.\r
 516       // -> See below paragraph 4.B, for the proposed treatment of such a case.\r
 517       \r
 518       // Add the vV0FromCascade to the aodVOVtxRefs\r
 519       fAODV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);\r
 520       \r
 521       \r
 522       // 4.A.2 - Add the positive tracks from the V0\r
 523       \r
 524       esdCascadePos->GetPxPyPz(momPos);\r
 525       esdCascadePos->GetXYZ(pos);\r
 526       esdCascadePos->GetCovarianceXYZPxPyPz(covTr);\r
 527       esdCascadePos->GetESDpid(pid);\r
 528       \r
 529       \r
 530       if (!fUsedTrack[idxPosFromV0Dghter]) {\r
 531         fUsedTrack[idxPosFromV0Dghter] = kTRUE;\r
 532         \r
 533         UInt_t selectInfo = 0;\r
 534         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);\r
 535         if(fMChandler) fMChandler->SelectParticle(esdCascadePos->GetLabel());\r
 536         aodTrack = new(tracksArray[fNumberOfTracks++]) \r
 537         AliAODTrack(  esdCascadePos->GetID(),\r
 538                     esdCascadePos->GetLabel(), \r
 539                     momPos, \r
 540                     kTRUE,\r
 541                     pos,\r
 542                     kFALSE, // Why kFALSE for "isDCA" ? FIXME\r
 543                     covTr, \r
 544                     (Short_t)esdCascadePos->GetSign(),\r
 545                     esdCascadePos->GetITSClusterMap(), \r
 546                     pid,\r
 547                     vV0FromCascade,\r
 548                     kTRUE,  // usedForVtxFit = kFALSE ? FIXME\r
 549                     vtx->UsesTrack(esdCascadePos->GetID()),\r
 550                     AliAODTrack::kSecondary,\r
 551                     selectInfo);\r
 552         aodTrack->SetTPCClusterMap(esdCascadePos->GetTPCClusterMap());\r
 553         aodTrack->SetTPCSharedMap (esdCascadePos->GetTPCSharedMap());\r
 554         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadePos));\r
 555         aodTrack->SetTPCPointsF(esdCascadePos->GetTPCNclsF());\r
 556         fAODTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);\r
 557         \r
 558         if (esdCascadePos->GetSign() > 0) ++fNumberOfPositiveTracks;\r
 559         aodTrack->ConvertAliPIDtoAODPID();\r
 560         aodTrack->SetFlags(esdCascadePos->GetStatus());\r
 561         SetAODPID(esdCascadePos,aodTrack,detpid);\r
 562       }\r
 563       else {\r
 564         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxPosFromV0Dghter));\r
 565       }\r
 566       vV0FromCascade->AddDaughter(aodTrack);\r
 567       \r
 568       \r
 569       // 4.A.3 - Add the negative tracks from the V0\r
 570       \r
 571       esdCascadeNeg->GetPxPyPz(momNeg);\r
 572       esdCascadeNeg->GetXYZ(pos);\r
 573       esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);\r
 574       esdCascadeNeg->GetESDpid(pid);\r
 575       \r
 576       \r
 577       if (!fUsedTrack[idxNegFromV0Dghter]) {\r
 578         fUsedTrack[idxNegFromV0Dghter] = kTRUE;\r
 579         \r
 580         UInt_t selectInfo = 0;\r
 581         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);\r
 582         if(fMChandler)fMChandler->SelectParticle(esdCascadeNeg->GetLabel());\r
 583         aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(  esdCascadeNeg->GetID(),\r
 584                                                       esdCascadeNeg->GetLabel(),\r
 585                                                       momNeg,\r
 586                                                       kTRUE,\r
 587                                                       pos,\r
 588                                                       kFALSE, // Why kFALSE for "isDCA" ? FIXME\r
 589                                                       covTr, \r
 590                                                       (Short_t)esdCascadeNeg->GetSign(),\r
 591                                                       esdCascadeNeg->GetITSClusterMap(), \r
 592                                                       pid,\r
 593                                                       vV0FromCascade,\r
 594                                                       kTRUE,  // usedForVtxFit = kFALSE ? FIXME\r
 595                                                       vtx->UsesTrack(esdCascadeNeg->GetID()),\r
 596                                                       AliAODTrack::kSecondary,\r
 597                                                       selectInfo);\r
 598         aodTrack->SetTPCClusterMap(esdCascadeNeg->GetTPCClusterMap());\r
 599         aodTrack->SetTPCSharedMap (esdCascadeNeg->GetTPCSharedMap());\r
 600         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeNeg));\r
 601         aodTrack->SetTPCPointsF(esdCascadeNeg->GetTPCNclsF());\r
 602         fAODTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);\r
 603         \r
 604         if (esdCascadeNeg->GetSign() > 0) ++fNumberOfPositiveTracks;\r
 605         aodTrack->ConvertAliPIDtoAODPID();\r
 606         aodTrack->SetFlags(esdCascadeNeg->GetStatus());\r
 607         SetAODPID(esdCascadeNeg,aodTrack,detpid);\r
 608       }\r
 609       else {\r
 610         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxNegFromV0Dghter));\r
 611       }\r
 612       \r
 613       vV0FromCascade->AddDaughter(aodTrack);\r
 614       \r
 615                         \r
 616       // 4.A.4 - Add the V0 from cascade to the V0 array\r
 617       \r
 618       Double_t  dcaV0Daughters      = esdV0FromCascade->GetDcaV0Daughters();\r
 619       Double_t  dcaV0ToPrimVertex   = esdV0FromCascade->GetD( esd.GetPrimaryVertex()->GetX(),\r
 620                                                              esd.GetPrimaryVertex()->GetY(),\r
 621                                                              esd.GetPrimaryVertex()->GetZ() );\r
 622       esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] ); \r
 623       esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] ); \r
 624       \r
 625       Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg\r
 626       dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD(      esd.GetPrimaryVertex()->GetX(),\r
 627                                                                   esd.GetPrimaryVertex()->GetY(),\r
 628                                                                   esd.GetMagneticField())        );\r
 629       dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD(      esd.GetPrimaryVertex()->GetX(),\r
 630                                                                   esd.GetPrimaryVertex()->GetY(),\r
 631                                                                   esd.GetMagneticField())        );\r
 632       \r
 633       aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0( vV0FromCascade, \r
 634                                                   dcaV0Daughters,\r
 635                                                   dcaV0ToPrimVertex, \r
 636                                                   momPosAtV0vtx, \r
 637                                                   momNegAtV0vtx, \r
 638                                                   dcaDaughterToPrimVertex); \r
 639       // set the aod v0 on-the-fly status\r
 640       aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());\r
 641       \r
 642       // Add the aodV0 to the aodVORefs\r
 643       fAODV0Refs->AddAt(aodV0,idxV0FromCascade);\r
 644       \r
 645       fUsedV0[idxV0FromCascade] = kTRUE;\r
 646       \r
 647     } else { \r
 648       // 4.B - if V0 structure already used\r
 649       \r
 650       // Note :\r
 651       //    one V0 can be used by several cascades (frequent in PbPb evts) : \r
 652       // same V0 which used but attached to different bachelor tracks\r
 653       // -> aodVORefs and fAODV0VtxRefs are needed.\r
 654       // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.\r
 655       \r
 656       vV0FromCascade = static_cast<AliAODVertex*>( fAODV0VtxRefs->At(idxV0FromCascade) );\r
 657       aodV0          = static_cast<AliAODv0*>    ( fAODV0Refs   ->At(idxV0FromCascade) );\r
 658       \r
 659       // - Treatment of the parent for such a "re-used" V0 :\r
 660       // Insert the cascade that reuses the V0 vertex in the lineage chain\r
 661       // Before : vV0 -> vCascade1 -> vPrimary\r
 662       //  - Hyp : cascade2 uses the same V0 as cascade1\r
 663       //  After :  vV0 -> vCascade2 -> vCascade1 -> vPrimary\r
 664       \r
 665       AliAODVertex *vCascadePreviousParent = static_cast<AliAODVertex*> (vV0FromCascade->GetParent());\r
 666                         vV0FromCascade->SetParent(vCascade);\r
 667                         vCascade      ->SetParent(vCascadePreviousParent);\r
 668       \r
 669 //      if(fDebug > 2)  \r
 670 //        printf("---- Cascade / Lineage insertion\n"\r
 671 //               "Parent of V0 vtx                 = Cascade vtx %p\n"\r
 672 //               "Parent of the cascade vtx        = Cascade vtx %p\n"\r
 673 //               "Parent of the parent cascade vtx = Cascade vtx %p\n", \r
 674 //               static_cast<void*> (vV0FromCascade->GetParent()),\r
 675 //               static_cast<void*> (vCascade->GetParent()),\r
 676 //               static_cast<void*> (vCascadePreviousParent->GetParent()) );\r
 677       \r
 678     }// end if V0 structure already used\r
 679     \r
 680 //    if (fDebug > 2) {\r
 681 //      printf("---- Cascade / V0 vertex: \n");\r
 682 //      vV0FromCascade->Print();\r
 683 //    }\r
 684 //    \r
 685 //    if (fDebug > 4) {\r
 686 //      printf("---- Cascade / pos dghter : \n");\r
 687 //                      aodTrack->Print();\r
 688 //      printf("---- Cascade / neg dghter : \n");\r
 689 //                      aodTrack->Print();\r
 690 //      printf("---- Cascade / aodV0 : \n");\r
 691 //                      aodV0->Print();\r
 692 //    }\r
 693     \r
 694     // In any case (used V0 or not), add the V0 vertex to the cascade one.\r
 695     vCascade->AddDaughter(vV0FromCascade);      \r
 696     \r
 697                 \r
 698     // 5 - Add the primary track of the cascade (if any)\r
 699     \r
 700     \r
 701     // 6 - Add the cascade to the AOD array of cascades\r
 702     \r
 703     Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd.GetPrimaryVertex()->GetX(),\r
 704                                                                      esd.GetPrimaryVertex()->GetY(),\r
 705                                                                      esd.GetMagneticField())        );\r
 706     \r
 707     Double_t momBachAtCascadeVtx[3]={0.};\r
 708 \r
 709     esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);\r
 710     \r
 711     aodCascade = new(cascadesArray[fNumberOfCascades++]) AliAODcascade(  vCascade,\r
 712                                                           esdCascade->Charge(),\r
 713                                                           esdCascade->GetDcaXiDaughters(),\r
 714                                                           -999.,\r
 715                                                           // DCAXiToPrimVtx -> needs to be calculated   ----|\r
 716                                                           // doesn't exist at ESD level;\r
 717                                                           // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)\r
 718                                                           dcaBachToPrimVertexXY,\r
 719                                                           momBachAtCascadeVtx,\r
 720                                                           *aodV0);\r
 721     \r
 722     if (fDebug > 10) {\r
 723       printf("---- Cascade / AOD cascade : \n\n");\r
 724       aodCascade->PrintXi(fPrimaryVertex->GetX(), fPrimaryVertex->GetY(), fPrimaryVertex->GetZ());\r
 725     }\r
 726     \r
 727   } // end of the loop on cascades\r
 728   \r
 729   Cascades().Expand(fNumberOfCascades);\r
 730 }\r
 731 \r
 732 //______________________________________________________________________________\r
 733 void AliAnalysisTaskESDfilter::ConvertV0s(const AliESDEvent& esd)\r
 734 {\r
 735   // Access to the AOD container of V0s\r
 736   \r
 737   AliCodeTimerAuto("",0);\r
 738 \r
 739   //\r
 740   // V0s\r
 741   //\r
 742   \r
 743   Double_t pos[3] = { 0. };      \r
 744   Double_t chi2(0.0);\r
 745   Double_t covVtx[6] = { 0. };\r
 746   Double_t momPos[3]={0.};\r
 747   Double_t covTr[21]={0.};\r
 748   Double_t pid[10]={0.};\r
 749   AliAODTrack* aodTrack(0x0);\r
 750   AliAODPid* detpid(0x0);\r
 751   Double_t momNeg[3]={0.};\r
 752   Double_t momPosAtV0vtx[3]={0.};\r
 753   Double_t momNegAtV0vtx[3]={0.};\r
 754     \r
 755   for (Int_t nV0 = 0; nV0 < esd.GetNumberOfV0s(); ++nV0) \r
 756   {\r
 757     if (fUsedV0[nV0]) continue; // skip if already added to the AOD\r
 758     \r
 759     AliESDv0 *v0 = esd.GetV0(nV0);\r
 760     Int_t posFromV0 = v0->GetPindex();\r
 761     Int_t negFromV0 = v0->GetNindex();\r
 762     \r
 763     // V0 selection \r
 764     //\r
 765     AliESDVertex *esdVtx   = new AliESDVertex(*(esd.GetPrimaryVertex()));\r
 766     AliESDtrack  *esdV0Pos = esd.GetTrack(posFromV0);\r
 767     AliESDtrack  *esdV0Neg = esd.GetTrack(negFromV0);\r
 768     TList v0objects;\r
 769     v0objects.AddAt(v0,                      0);\r
 770     v0objects.AddAt(esdV0Pos,                1);\r
 771     v0objects.AddAt(esdV0Neg,                2);\r
 772     v0objects.AddAt(esdVtx,                  3);\r
 773     UInt_t selectV0 = 0;\r
 774     if (fV0Filter) {\r
 775       selectV0 = fV0Filter->IsSelected(&v0objects);\r
 776       // this is a little awkward but otherwise the \r
 777       // list wants to access the pointer (delete it) \r
 778       // again when going out of scope\r
 779       delete v0objects.RemoveAt(3); // esdVtx created via copy construct\r
 780       esdVtx = 0;\r
 781       if (!selectV0) \r
 782         continue;\r
 783     }\r
 784     else{\r
 785       delete v0objects.RemoveAt(3); // esdVtx created via copy construct\r
 786       esdVtx = 0;\r
 787     }\r
 788     \r
 789     v0->GetXYZ(pos[0], pos[1], pos[2]);\r
 790     \r
 791     if (!fOldESDformat) {\r
 792             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3\r
 793             v0->GetPosCov(covVtx);\r
 794     } else {\r
 795             chi2 = -999.;\r
 796             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;\r
 797     }\r
 798     \r
 799     \r
 800     AliAODVertex * vV0 = \r
 801           new(Vertices()[fNumberOfVertices++]) AliAODVertex(pos,\r
 802                                                       covVtx,\r
 803                                                       chi2,\r
 804                                                       fPrimaryVertex,\r
 805                                                       nV0,\r
 806                                                       AliAODVertex::kV0);\r
 807     fPrimaryVertex->AddDaughter(vV0);\r
 808     \r
 809     \r
 810     // Add the positive tracks from the V0\r
 811     \r
 812 \r
 813     esdV0Pos->GetPxPyPz(momPos);\r
 814     esdV0Pos->GetXYZ(pos);\r
 815     esdV0Pos->GetCovarianceXYZPxPyPz(covTr);\r
 816     esdV0Pos->GetESDpid(pid);\r
 817     \r
 818     const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
 819 \r
 820     if (!fUsedTrack[posFromV0]) {\r
 821             fUsedTrack[posFromV0] = kTRUE;\r
 822             UInt_t selectInfo = 0;\r
 823             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);\r
 824             if(fMChandler)fMChandler->SelectParticle(esdV0Pos->GetLabel());\r
 825             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Pos->GetID(),\r
 826                                                     esdV0Pos->GetLabel(), \r
 827                                                     momPos, \r
 828                                                     kTRUE,\r
 829                                                     pos,\r
 830                                                     kFALSE,\r
 831                                                     covTr, \r
 832                                                     (Short_t)esdV0Pos->GetSign(),\r
 833                                                     esdV0Pos->GetITSClusterMap(), \r
 834                                                     pid,\r
 835                                                     vV0,\r
 836                                                     kTRUE,  // check if this is right\r
 837                                                     vtx->UsesTrack(esdV0Pos->GetID()),\r
 838                                                     AliAODTrack::kSecondary,\r
 839                                                     selectInfo);\r
 840             aodTrack->SetTPCClusterMap(esdV0Pos->GetTPCClusterMap());\r
 841             aodTrack->SetTPCSharedMap (esdV0Pos->GetTPCSharedMap());\r
 842             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Pos));\r
 843             aodTrack->SetTPCPointsF(esdV0Pos->GetTPCNclsF());\r
 844             fAODTrackRefs->AddAt(aodTrack,posFromV0);\r
 845             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());\r
 846             if (esdV0Pos->GetSign() > 0) ++fNumberOfPositiveTracks;\r
 847             aodTrack->ConvertAliPIDtoAODPID();\r
 848             aodTrack->SetFlags(esdV0Pos->GetStatus());\r
 849       SetAODPID(esdV0Pos,aodTrack,detpid);\r
 850     }\r
 851     else {\r
 852             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));\r
 853             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());\r
 854     }\r
 855     vV0->AddDaughter(aodTrack);\r
 856     \r
 857     // Add the negative tracks from the V0\r
 858     \r
 859     esdV0Neg->GetPxPyPz(momNeg);\r
 860     esdV0Neg->GetXYZ(pos);\r
 861     esdV0Neg->GetCovarianceXYZPxPyPz(covTr);\r
 862     esdV0Neg->GetESDpid(pid);\r
 863     \r
 864     if (!fUsedTrack[negFromV0]) {\r
 865             fUsedTrack[negFromV0] = kTRUE;\r
 866             UInt_t selectInfo = 0;\r
 867             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);\r
 868             if(fMChandler)fMChandler->SelectParticle(esdV0Neg->GetLabel());\r
 869             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Neg->GetID(),\r
 870                                                     esdV0Neg->GetLabel(),\r
 871                                                     momNeg,\r
 872                                                     kTRUE,\r
 873                                                     pos,\r
 874                                                     kFALSE,\r
 875                                                     covTr, \r
 876                                                     (Short_t)esdV0Neg->GetSign(),\r
 877                                                     esdV0Neg->GetITSClusterMap(), \r
 878                                                     pid,\r
 879                                                     vV0,\r
 880                                                     kTRUE,  // check if this is right\r
 881                                                     vtx->UsesTrack(esdV0Neg->GetID()),\r
 882                                                     AliAODTrack::kSecondary,\r
 883                                                     selectInfo);\r
 884             aodTrack->SetTPCClusterMap(esdV0Neg->GetTPCClusterMap());\r
 885             aodTrack->SetTPCSharedMap (esdV0Neg->GetTPCSharedMap());\r
 886             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Neg));\r
 887             aodTrack->SetTPCPointsF(esdV0Neg->GetTPCNclsF());\r
 888             \r
 889             fAODTrackRefs->AddAt(aodTrack,negFromV0);\r
 890             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());\r
 891             if (esdV0Neg->GetSign() > 0) ++fNumberOfPositiveTracks;\r
 892             aodTrack->ConvertAliPIDtoAODPID();\r
 893             aodTrack->SetFlags(esdV0Neg->GetStatus());\r
 894       SetAODPID(esdV0Neg,aodTrack,detpid);\r
 895     }\r
 896     else {\r
 897             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(negFromV0));\r
 898             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());\r
 899     }\r
 900     vV0->AddDaughter(aodTrack);\r
 901     \r
 902     \r
 903     // Add the V0 the V0 array as well\r
 904     \r
 905     Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();\r
 906     Double_t  dcaV0ToPrimVertex   = v0->GetD(esd.GetPrimaryVertex()->GetX(),\r
 907                                              esd.GetPrimaryVertex()->GetY(),\r
 908                                              esd.GetPrimaryVertex()->GetZ());\r
 909     \r
 910     v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]); \r
 911     v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]); \r
 912     \r
 913     Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg\r
 914     dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD(  esd.GetPrimaryVertex()->GetX(),\r
 915                                                            esd.GetPrimaryVertex()->GetY(),\r
 916                                                            esd.GetMagneticField()) );\r
 917     dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD(  esd.GetPrimaryVertex()->GetX(),\r
 918                                                            esd.GetPrimaryVertex()->GetY(),\r
 919                                                            esd.GetMagneticField()) );\r
 920     \r
 921     AliAODv0* aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0(vV0, \r
 922                                                 dcaV0Daughters,\r
 923                                                 dcaV0ToPrimVertex,\r
 924                                                 momPosAtV0vtx,\r
 925                                                 momNegAtV0vtx,\r
 926                                                 dcaDaughterToPrimVertex);\r
 927     \r
 928     // set the aod v0 on-the-fly status\r
 929     aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());\r
 930   }//End of loop on V0s \r
 931   \r
 932   V0s().Expand(fNumberOfV0s);    \r
 933 }\r
 934 \r
 935 //______________________________________________________________________________\r
 936 void AliAnalysisTaskESDfilter::ConvertTPCOnlyTracks(const AliESDEvent& esd)\r
 937 {\r
 938   // Convert TPC only tracks\r
 939   \r
 940   AliCodeTimerAuto("",0);\r
 941   \r
 942   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks\r
 943   for(int it = 0;it < fNumberOfTracks;++it)\r
 944   {\r
 945     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));\r
 946     if(!tr)continue;\r
 947     UInt_t map = tr->GetFilterMap();\r
 948     if(map&fTPCOnlyFilterMask){\r
 949       // we only reset the track select ionfo, no deletion...\r
 950       tr->SetFilterMap(map&~fTPCOnlyFilterMask);\r
 951     }\r
 952   }\r
 953   // Loop over the ESD trcks and pick out the tracks passing TPC only cuts\r
 954   \r
 955   \r
 956   const AliESDVertex *vtxSPD = esd.GetPrimaryVertexSPD();\r
 957   const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
 958 \r
 959   Double_t pos[3] = { 0. };      \r
 960   Double_t covTr[21]={0.};\r
 961   Double_t pid[10]={0.};  \r
 962   Double_t p[3] = { 0. };\r
 963 \r
 964   Double_t pDCA[3] = { 0. }; // momentum at DCA\r
 965   Double_t rDCA[3] = { 0. }; // position at DCA\r
 966   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z\r
 967   Float_t  cDCA[3] = {0.};    // covariance of impact parameters\r
 968 \r
 969 \r
 970   AliAODTrack* aodTrack(0x0);\r
 971   \r
 972   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) \r
 973   {\r
 974     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy \r
 975     \r
 976     UInt_t selectInfo = 0;\r
 977     //\r
 978     // Track selection\r
 979     if (fTrackFilter) {\r
 980       selectInfo = fTrackFilter->IsSelected(esdTrack);\r
 981     }\r
 982     selectInfo &= fTPCOnlyFilterMask;\r
 983     if (!selectInfo)continue;\r
 984     \r
 985     // create a tpc only tracl\r
 986     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(&esd),esdTrack->GetID());\r
 987     if(!track) continue;\r
 988     \r
 989     if(track->Pt()>0.)\r
 990     {\r
 991       // only constrain tracks above threshold\r
 992       AliExternalTrackParam exParam;\r
 993       // take the B-field from the ESD, no 3D fieldMap available at this point\r
 994       Bool_t relate = false;\r
 995       relate = track->RelateToVertexTPC(vtxSPD,esd.GetMagneticField(),kVeryBig,&exParam);\r
 996       if(!relate){\r
 997         delete track;\r
 998         continue;\r
 999       }\r
1000       // fetch the track parameters at the DCA (unconstraint)\r
1001       if(track->GetTPCInnerParam()){\r
1002         track->GetTPCInnerParam()->GetPxPyPz(pDCA);\r
1003         track->GetTPCInnerParam()->GetXYZ(rDCA);\r
1004       }\r
1005       // get the DCA to the vertex:\r
1006       track->GetImpactParametersTPC(dDCA,cDCA);\r
1007       // set the constraint parameters to the track\r
1008       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());\r
1009     }\r
1010     \r
1011     track->GetPxPyPz(p);\r
1012     track->GetXYZ(pos);\r
1013     track->GetCovarianceXYZPxPyPz(covTr);\r
1014     track->GetESDpid(pid);\r
1015     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());\r
1016     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((track->GetID()+1)*-1,\r
1017                                                             track->GetLabel(),\r
1018                                                             p,\r
1019                                                             kTRUE,\r
1020                                                             pos,\r
1021                                                             kFALSE,\r
1022                                                             covTr, \r
1023                                                             (Short_t)track->GetSign(),\r
1024                                                             track->GetITSClusterMap(), \r
1025                                                             pid,\r
1026                                                             fPrimaryVertex,\r
1027                                                             kTRUE, // check if this is right\r
1028                                                             vtx->UsesTrack(track->GetID()),\r
1029                                                             AliAODTrack::kPrimary, \r
1030                                                             selectInfo);\r
1031     aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());\r
1032     aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());\r
1033     Float_t ndf = track->GetTPCNcls()+1 - 5 ;\r
1034     if(ndf>0){\r
1035       aodTrack->SetChi2perNDF(track->GetConstrainedChi2TPC());\r
1036     }\r
1037     else{\r
1038       aodTrack->SetChi2perNDF(-1);\r
1039     }\r
1040 \r
1041     // set the DCA values to the AOD track\r
1042     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);\r
1043     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);\r
1044     aodTrack->SetDCA(dDCA[0],dDCA[1]);\r
1045 \r
1046     aodTrack->SetFlags(track->GetStatus());\r
1047     aodTrack->SetTPCPointsF(track->GetTPCNclsF());\r
1048 \r
1049     delete track;\r
1050   } // end of loop on tracks\r
1051   \r
1052 }\r
1053 \r
1054 //______________________________________________________________________________\r
1055 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)\r
1056 {\r
1057   // Tracks (primary and orphan)\r
1058 \r
1059   AliCodeTimerAuto("",0);\r
1060   \r
1061   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));\r
1062   \r
1063   const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
1064   Double_t p[3] = { 0. };\r
1065   Double_t pos[3] = { 0. };\r
1066   Double_t covTr[21] = { 0. };\r
1067   Double_t pid[10] = { 0. };\r
1068   AliAODTrack* aodTrack(0x0);\r
1069   AliAODPid* detpid(0x0);\r
1070   \r
1071   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) \r
1072   {\r
1073     if (fUsedTrack[nTrack]) continue;\r
1074     \r
1075     AliESDtrack *esdTrack = esd.GetTrack(nTrack);\r
1076     UInt_t selectInfo = 0;\r
1077     //\r
1078     // Track selection\r
1079     if (fTrackFilter) {\r
1080             selectInfo = fTrackFilter->IsSelected(esdTrack);\r
1081             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;\r
1082     }\r
1083     \r
1084     \r
1085     esdTrack->GetPxPyPz(p);\r
1086     esdTrack->GetXYZ(pos);\r
1087     esdTrack->GetCovarianceXYZPxPyPz(covTr);\r
1088     esdTrack->GetESDpid(pid);\r
1089     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());\r
1090     fPrimaryVertex->AddDaughter(aodTrack =\r
1091                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),\r
1092                                                             esdTrack->GetLabel(),\r
1093                                                             p,\r
1094                                                             kTRUE,\r
1095                                                             pos,\r
1096                                                             kFALSE,\r
1097                                                             covTr, \r
1098                                                             (Short_t)esdTrack->GetSign(),\r
1099                                                             esdTrack->GetITSClusterMap(), \r
1100                                                             pid,\r
1101                                                             fPrimaryVertex,\r
1102                                                             kTRUE, // check if this is right\r
1103                                                             vtx->UsesTrack(esdTrack->GetID()),\r
1104                                                             AliAODTrack::kPrimary, \r
1105                                                             selectInfo)\r
1106                          );\r
1107     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());\r
1108     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());\r
1109     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));\r
1110     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());\r
1111 \r
1112     fAODTrackRefs->AddAt(aodTrack, nTrack);\r
1113     \r
1114     \r
1115     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1116     aodTrack->SetFlags(esdTrack->GetStatus());\r
1117     aodTrack->ConvertAliPIDtoAODPID();\r
1118     SetAODPID(esdTrack,aodTrack,detpid);\r
1119   } // end of loop on tracks\r
1120 }\r
1121 \r
1122 //______________________________________________________________________________\r
1123 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)\r
1124 {\r
1125   AliCodeTimerAuto("",0);\r
1126   Int_t jPmdClusters=0;\r
1127   // Access to the AOD container of PMD clusters\r
1128   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());\r
1129   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {\r
1130     // file pmd clusters, to be revised!\r
1131     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);\r
1132     Int_t nLabel = 0;\r
1133     Int_t *label = 0x0;\r
1134     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};\r
1135     Double_t pidPmd[13] = { 0.}; // to be revised!\r
1136     // type not set!\r
1137     // assoc cluster not set\r
1138     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);\r
1139   }\r
1140 }\r
1141 \r
1142 //______________________________________________________________________________\r
1143 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)\r
1144 {\r
1145   AliCodeTimerAuto("",0);\r
1146   \r
1147   // Access to the AOD container of clusters\r
1148   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());\r
1149   Int_t jClusters(0);\r
1150   \r
1151   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {\r
1152     \r
1153     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);\r
1154     \r
1155     Int_t  id        = cluster->GetID();\r
1156     Int_t  nLabel    = cluster->GetNLabels();\r
1157     Int_t *labels    = cluster->GetLabels();\r
1158     if(labels){ \r
1159                   for(int i = 0;i < nLabel;++i){\r
1160                           if(fMChandler)fMChandler->SelectParticle(labels[i]);\r
1161                   }\r
1162           }             \r
1163     \r
1164     Float_t energy = cluster->E();\r
1165     Float_t posF[3] = { 0.};\r
1166     cluster->GetPosition(posF);\r
1167     \r
1168     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,\r
1169                                                                                       nLabel,\r
1170                                                                                       labels,\r
1171                                                                                       energy,\r
1172                                                                                       posF,\r
1173                                                                                       NULL,\r
1174                                                                                       cluster->GetType(),0);\r
1175     \r
1176     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),\r
1177                                 cluster->GetDispersion(),\r
1178                                 cluster->GetM20(), cluster->GetM02(),\r
1179                                 cluster->GetEmcCpvDistance(),  \r
1180                                 cluster->GetNExMax(),cluster->GetTOF()) ;\r
1181     \r
1182     caloCluster->SetPIDFromESD(cluster->GetPID());\r
1183     caloCluster->SetNCells(cluster->GetNCells());\r
1184     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());\r
1185     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());\r
1186     \r
1187     TArrayI* matchedT =         cluster->GetTracksMatched();\r
1188     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {        \r
1189       for (Int_t im = 0; im < matchedT->GetSize(); im++) {\r
1190         Int_t iESDtrack = matchedT->At(im);;\r
1191         if (fAODTrackRefs->At(iESDtrack) != 0) {\r
1192           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));\r
1193         }\r
1194       }\r
1195     }\r
1196     \r
1197   } \r
1198   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      \r
1199 }\r
1200 \r
1201 //______________________________________________________________________________\r
1202 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)\r
1203 {\r
1204   AliCodeTimerAuto("",0);\r
1205   // fill EMCAL cell info\r
1206   if (esd.GetEMCALCells()) { // protection against missing ESD information\r
1207     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());\r
1208     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;\r
1209     \r
1210     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());\r
1211     aodEMcells.CreateContainer(nEMcell);\r
1212     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);\r
1213     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      \r
1214       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell));\r
1215     }\r
1216     aodEMcells.Sort();\r
1217   }\r
1218 }\r
1219 \r
1220 //______________________________________________________________________________\r
1221 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)\r
1222 {\r
1223   AliCodeTimerAuto("",0);\r
1224   // fill PHOS cell info\r
1225   if (esd.GetPHOSCells()) { // protection against missing ESD information\r
1226     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());\r
1227     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;\r
1228     \r
1229     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());\r
1230     aodPHcells.CreateContainer(nPHcell);\r
1231     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);\r
1232     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      \r
1233       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));\r
1234     }\r
1235     aodPHcells.Sort();\r
1236   }\r
1237 }\r
1238 \r
1239 //______________________________________________________________________________\r
1240 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)\r
1241 {\r
1242   // tracklets    \r
1243   AliCodeTimerAuto("",0);\r
1244 \r
1245   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());\r
1246   const AliMultiplicity *mult = esd.GetMultiplicity();\r
1247   if (mult) {\r
1248     if (mult->GetNumberOfTracklets()>0) {\r
1249       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());\r
1250       \r
1251       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {\r
1252         if(fMChandler){\r
1253           fMChandler->SelectParticle(mult->GetLabel(n, 0));\r
1254           fMChandler->SelectParticle(mult->GetLabel(n, 1));\r
1255         }\r
1256         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));\r
1257       }\r
1258     }\r
1259   } else {\r
1260     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");\r
1261   }\r
1262 }\r
1263 \r
1264 //______________________________________________________________________________\r
1265 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)\r
1266 {\r
1267   AliCodeTimerAuto("",0);\r
1268   \r
1269   // Kinks: it is a big mess the access to the information in the kinks\r
1270   // The loop is on the tracks in order to find the mother and daugther of each kink\r
1271   \r
1272   Double_t covTr[21]={0.};\r
1273   Double_t pid[10]={0.};\r
1274   AliAODPid* detpid(0x0);\r
1275   \r
1276   fNumberOfKinks = esd.GetNumberOfKinks();\r
1277 \r
1278   const AliESDVertex* vtx = esd.GetPrimaryVertex();\r
1279   \r
1280   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack) \r
1281   {\r
1282     AliESDtrack * esdTrack = esd.GetTrack(iTrack);\r
1283     \r
1284     Int_t ikink = esdTrack->GetKinkIndex(0);\r
1285     \r
1286     if (ikink && fNumberOfKinks) {\r
1287             // Negative kink index: mother, positive: daughter\r
1288             \r
1289             // Search for the second track of the kink\r
1290             \r
1291             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {\r
1292         \r
1293         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);\r
1294         \r
1295         Int_t jkink = esdTrack1->GetKinkIndex(0);\r
1296         \r
1297         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {\r
1298           \r
1299           // The two tracks are from the same kink\r
1300           \r
1301           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks\r
1302           \r
1303           Int_t imother = -1;\r
1304           Int_t idaughter = -1;\r
1305           \r
1306           if (ikink<0 && jkink>0) {\r
1307             \r
1308             imother = iTrack;\r
1309             idaughter = jTrack;\r
1310           }\r
1311           else if (ikink>0 && jkink<0) {\r
1312             \r
1313             imother = jTrack;\r
1314             idaughter = iTrack;\r
1315           }\r
1316           else {\r
1317             //                  cerr << "Error: Wrong combination of kink indexes: "\r
1318             //                       << ikink << " " << jkink << endl;\r
1319             continue;\r
1320           }\r
1321           \r
1322           // Add the mother track if it passed primary track selection cuts\r
1323           \r
1324           AliAODTrack * mother = NULL;\r
1325           \r
1326           UInt_t selectInfo = 0;\r
1327           if (fTrackFilter) {\r
1328             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));\r
1329             if (!selectInfo) continue;\r
1330           }\r
1331           \r
1332           if (!fUsedTrack[imother]) {\r
1333             \r
1334             fUsedTrack[imother] = kTRUE;\r
1335             \r
1336             AliESDtrack *esdTrackM = esd.GetTrack(imother);\r
1337             Double_t p[3] = { 0. };\r
1338             Double_t pos[3] = { 0. };\r
1339             esdTrackM->GetPxPyPz(p);\r
1340             esdTrackM->GetXYZ(pos);\r
1341             esdTrackM->GetCovarianceXYZPxPyPz(covTr);\r
1342             esdTrackM->GetESDpid(pid);\r
1343             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());\r
1344             mother = \r
1345             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),\r
1346                                                esdTrackM->GetLabel(),\r
1347                                                p,\r
1348                                                kTRUE,\r
1349                                                pos,\r
1350                                                kFALSE,\r
1351                                                covTr, \r
1352                                                (Short_t)esdTrackM->GetSign(),\r
1353                                                esdTrackM->GetITSClusterMap(), \r
1354                                                pid,\r
1355                                                fPrimaryVertex,\r
1356                                                kTRUE, // check if this is right\r
1357                                                vtx->UsesTrack(esdTrack->GetID()),\r
1358                                                AliAODTrack::kPrimary,\r
1359                                                selectInfo);\r
1360             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());\r
1361             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());\r
1362             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));\r
1363             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());\r
1364 \r
1365             fAODTrackRefs->AddAt(mother, imother);\r
1366             \r
1367             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1368             mother->SetFlags(esdTrackM->GetStatus());\r
1369             mother->ConvertAliPIDtoAODPID();\r
1370             fPrimaryVertex->AddDaughter(mother);\r
1371             mother->ConvertAliPIDtoAODPID();\r
1372             SetAODPID(esdTrackM,mother,detpid);\r
1373           }\r
1374           else {\r
1375             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1\r
1376             //                       << " track " << imother << " has already been used!" << endl;\r
1377           }\r
1378           \r
1379           // Add the kink vertex\r
1380           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);\r
1381           \r
1382           AliAODVertex * vkink = \r
1383           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),\r
1384                                                   NULL,\r
1385                                                   0.,\r
1386                                                   mother,\r
1387                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!\r
1388                                                   AliAODVertex::kKink);\r
1389           // Add the daughter track\r
1390           \r
1391           AliAODTrack * daughter = NULL;\r
1392           \r
1393           if (!fUsedTrack[idaughter]) {\r
1394             \r
1395             fUsedTrack[idaughter] = kTRUE;\r
1396             \r
1397             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);\r
1398             Double_t p[3] = { 0. };\r
1399             Double_t pos[3] = { 0. };\r
1400 \r
1401             esdTrackD->GetPxPyPz(p);\r
1402             esdTrackD->GetXYZ(pos);\r
1403             esdTrackD->GetCovarianceXYZPxPyPz(covTr);\r
1404             esdTrackD->GetESDpid(pid);\r
1405             selectInfo = 0;\r
1406             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);\r
1407             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());\r
1408             daughter = \r
1409             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),\r
1410                                                esdTrackD->GetLabel(),\r
1411                                                p,\r
1412                                                kTRUE,\r
1413                                                pos,\r
1414                                                kFALSE,\r
1415                                                covTr, \r
1416                                                (Short_t)esdTrackD->GetSign(),\r
1417                                                esdTrackD->GetITSClusterMap(), \r
1418                                                pid,\r
1419                                                vkink,\r
1420                                                kTRUE, // check if this is right\r
1421                                                vtx->UsesTrack(esdTrack->GetID()),\r
1422                                                AliAODTrack::kSecondary,\r
1423                                                selectInfo);\r
1424             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());\r
1425             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());\r
1426             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());\r
1427             fAODTrackRefs->AddAt(daughter, idaughter);\r
1428             \r
1429             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1430             daughter->SetFlags(esdTrackD->GetStatus());\r
1431             daughter->ConvertAliPIDtoAODPID();\r
1432             vkink->AddDaughter(daughter);\r
1433             daughter->ConvertAliPIDtoAODPID();\r
1434             SetAODPID(esdTrackD,daughter,detpid);\r
1435           }\r
1436           else {\r
1437             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1\r
1438             //                       << " track " << idaughter << " has already been used!" << endl;\r
1439           }\r
1440         }\r
1441             }\r
1442     }      \r
1443   }\r
1444 }\r
1445 \r
1446 //______________________________________________________________________________\r
1447 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)\r
1448 {\r
1449   AliCodeTimerAuto("",0);\r
1450   \r
1451   // Access to the AOD container of vertices\r
1452   fNumberOfVertices = 0;\r
1453   \r
1454   Double_t pos[3] = { 0. };\r
1455   Double_t covVtx[6] = { 0. };\r
1456 \r
1457   // Add primary vertex. The primary tracks will be defined\r
1458   // after the loops on the composite objects (V0, cascades, kinks)\r
1459   const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
1460   \r
1461   vtx->GetXYZ(pos); // position\r
1462   vtx->GetCovMatrix(covVtx); //covariance matrix\r
1463   \r
1464   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])\r
1465   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);\r
1466   fPrimaryVertex->SetName(vtx->GetName());\r
1467   fPrimaryVertex->SetTitle(vtx->GetTitle());\r
1468   \r
1469   TString vtitle = vtx->GetTitle();\r
1470   if (!vtitle.Contains("VertexerTracks")) \r
1471     fPrimaryVertex->SetNContributors(vtx->GetNContributors());\r
1472   \r
1473   if (fDebug > 0) fPrimaryVertex->Print();  \r
1474   \r
1475   // Add SPD "main" vertex \r
1476   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();\r
1477   vtxS->GetXYZ(pos); // position\r
1478   vtxS->GetCovMatrix(covVtx); //covariance matrix\r
1479   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])\r
1480   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);\r
1481   mVSPD->SetName(vtxS->GetName());\r
1482   mVSPD->SetTitle(vtxS->GetTitle());\r
1483   mVSPD->SetNContributors(vtxS->GetNContributors()); \r
1484   \r
1485   // Add SPD pileup vertices\r
1486   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)\r
1487   {\r
1488     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);\r
1489     vtxP->GetXYZ(pos); // position\r
1490     vtxP->GetCovMatrix(covVtx); //covariance matrix\r
1491     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])\r
1492     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);\r
1493     pVSPD->SetName(vtxP->GetName());\r
1494     pVSPD->SetTitle(vtxP->GetTitle());\r
1495     pVSPD->SetNContributors(vtxP->GetNContributors()); \r
1496   }\r
1497   \r
1498   // Add TRK pileup vertices\r
1499   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)\r
1500   {\r
1501     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);\r
1502     vtxP->GetXYZ(pos); // position\r
1503     vtxP->GetCovMatrix(covVtx); //covariance matrix\r
1504     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])\r
1505     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);\r
1506     pVTRK->SetName(vtxP->GetName());\r
1507     pVTRK->SetTitle(vtxP->GetTitle());\r
1508     pVTRK->SetNContributors(vtxP->GetNContributors());\r
1509   }\r
1510 }\r
1511 \r
1512 //______________________________________________________________________________\r
1513 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)\r
1514 {\r
1515   // Convert VZERO data\r
1516   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();\r
1517   *vzeroData = *(esd.GetVZEROData());\r
1518 }\r
1519 \r
1520 //______________________________________________________________________________\r
1521 void AliAnalysisTaskESDfilter::ConvertESDtoAOD() \r
1522 {\r
1523   // ESD Filter analysis task executed for each event\r
1524   \r
1525   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
1526   \r
1527   if(!esd)return;\r
1528 \r
1529   AliCodeTimerAuto("",0);\r
1530   \r
1531   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );\r
1532   \r
1533   fNumberOfTracks = 0;\r
1534   fNumberOfPositiveTracks = 0;\r
1535   fNumberOfV0s = 0;\r
1536   fNumberOfVertices = 0;\r
1537   fNumberOfCascades = 0;\r
1538   fNumberOfKinks = 0;\r
1539     \r
1540   AliAODHeader* header = ConvertHeader(*esd);\r
1541 \r
1542   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);\r
1543   \r
1544   // Fetch Stack for debuggging if available \r
1545   fMChandler=0x0;\r
1546   if(MCEvent())\r
1547   {\r
1548     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); \r
1549   }\r
1550   \r
1551   // loop over events and fill them\r
1552   // Multiplicity information needed by the header (to be revised!)\r
1553   Int_t nTracks    = esd->GetNumberOfTracks();\r
1554   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);\r
1555 \r
1556   // Update the header\r
1557 \r
1558   Int_t nV0s      = esd->GetNumberOfV0s();\r
1559   Int_t nCascades = esd->GetNumberOfCascades();\r
1560   Int_t nKinks    = esd->GetNumberOfKinks();\r
1561   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;\r
1562   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex\r
1563   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();\r
1564   nVertices+=nPileSPDVertices;\r
1565   nVertices+=nPileTrkVertices;\r
1566   Int_t nJets     = 0;\r
1567   Int_t nCaloClus = esd->GetNumberOfCaloClusters();\r
1568   Int_t nFmdClus  = 0;\r
1569   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();\r
1570     \r
1571   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));\r
1572        \r
1573   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);\r
1574 \r
1575   if (nV0s > 0) \r
1576   {\r
1577     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0\r
1578     fAODV0VtxRefs = new TRefArray(nV0s);\r
1579     // RefArray to store the mapping between esd V0 number and newly created AOD-V0\r
1580     fAODV0Refs = new TRefArray(nV0s); \r
1581     // Array to take into account the V0s already added to the AOD (V0 within cascades)\r
1582     fUsedV0 = new Bool_t[nV0s];\r
1583     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;\r
1584   }\r
1585   \r
1586   if (nTracks>0) \r
1587   {\r
1588     // RefArray to store the mapping between esd track number and newly created AOD-Track\r
1589     \r
1590     fAODTrackRefs = new TRefArray(nTracks);\r
1591 \r
1592     // Array to take into account the tracks already added to the AOD    \r
1593     fUsedTrack = new Bool_t[nTracks];\r
1594     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;\r
1595   }\r
1596   \r
1597   // Array to take into account the kinks already added to the AOD\r
1598   if (nKinks>0) \r
1599   {\r
1600     fUsedKink = new Bool_t[nKinks];\r
1601     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;\r
1602   }\r
1603     \r
1604   ConvertPrimaryVertices(*esd);\r
1605 \r
1606   //setting best TOF PID\r
1607   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
1608   if (esdH)\r
1609       fESDpid = esdH->GetESDpid();\r
1610 \r
1611   if (fIsPidOwner && fESDpid){\r
1612     delete fESDpid;\r
1613     fESDpid = 0;\r
1614   }\r
1615   if(!fESDpid)\r
1616   { //in case of no Tender attached \r
1617     fESDpid = new AliESDpid;\r
1618     fIsPidOwner = kTRUE;\r
1619   }\r
1620   \r
1621   if(!esd->GetTOFHeader())\r
1622   { //protection in case the pass2 LHC10b,c,d have been processed without tender. \r
1623     Float_t t0spread[10];\r
1624     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!! \r
1625     for (Int_t i=0; i<10; i++) t0spread[i] = (TMath::Sqrt(esd->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps\r
1626     fESDpid->GetTOFResponse().SetT0resolution(t0spread);\r
1627     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);\r
1628     \r
1629     fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);    \r
1630   }\r
1631   \r
1632   if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender. \r
1633   \r
1634   if ( fAreCascadesEnabled ) ConvertCascades(*esd);\r
1635 \r
1636   if ( fAreV0sEnabled ) ConvertV0s(*esd);\r
1637   \r
1638   if ( fAreKinksEnabled ) ConvertKinks(*esd);\r
1639   \r
1640   if ( fAreTracksEnabled ) ConvertTracks(*esd);\r
1641   \r
1642   // Update number of AOD tracks in header at the end of track loop (M.G.)\r
1643   header->SetRefMultiplicity(fNumberOfTracks);\r
1644   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);\r
1645   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);\r
1646 \r
1647   if ( fTPCOnlyFilterMask ) ConvertTPCOnlyTracks(*esd);\r
1648   \r
1649   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);\r
1650   \r
1651   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);\r
1652   \r
1653   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);\r
1654   \r
1655   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);\r
1656   \r
1657   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);\r
1658   \r
1659   delete fAODTrackRefs; fAODTrackRefs=0x0;\r
1660   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;\r
1661   delete fAODV0Refs; fAODV0Refs=0x0;\r
1662   \r
1663   delete[] fUsedTrack; fUsedTrack=0x0;\r
1664   delete[] fUsedV0; fUsedV0=0x0;\r
1665   delete[] fUsedKink; fUsedKink=0x0;\r
1666 \r
1667   if ( fIsPidOwner){\r
1668     delete fESDpid;\r
1669     fESDpid = 0x0;\r
1670   }\r
1671 \r
1672 \r
1673 }\r
1674 \r
1675 \r
1676 //______________________________________________________________________________\r
1677 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)\r
1678 {\r
1679   //\r
1680   // Setter for the raw PID detector signals\r
1681   //\r
1682 \r
1683   // Save PID object for candidate electrons\r
1684     Bool_t pidSave = kFALSE;\r
1685     if (fTrackFilter) {\r
1686         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");\r
1687         if (selectInfo)  pidSave = kTRUE;\r
1688     }\r
1689 \r
1690 \r
1691     // Tracks passing pt cut \r
1692     if(esdtrack->Pt()>fHighPthreshold) {\r
1693         pidSave = kTRUE;\r
1694     } else {\r
1695         if(fPtshape){\r
1696             if(esdtrack->Pt()> fPtshape->GetXmin()){\r
1697                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);\r
1698                 if(gRandom->Rndm(0)<1./y){\r
1699                     pidSave = kTRUE;\r
1700                 }//end rndm\r
1701             }//end if p < pmin\r
1702         }//end if p function\r
1703     }// end else\r
1704 \r
1705     if (pidSave) {\r
1706       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track\r
1707         detpid = new AliAODPid();\r
1708         SetDetectorRawSignals(detpid,esdtrack);\r
1709         aodtrack->SetDetPID(detpid);\r
1710       }\r
1711     }\r
1712 }\r
1713 \r
1714 //______________________________________________________________________________\r
1715 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)\r
1716 {\r
1717 //\r
1718 //assignment of the detector signals (AliXXXesdPID inspired)\r
1719 //\r
1720  if(!track){\r
1721  AliInfo("no ESD track found. .....exiting");\r
1722  return;\r
1723  }\r
1724  // TPC momentum\r
1725  const AliExternalTrackParam *in=track->GetInnerParam();\r
1726  if (in) {\r
1727    aodpid->SetTPCmomentum(in->GetP());\r
1728  }else{\r
1729    aodpid->SetTPCmomentum(-1.);\r
1730  }\r
1731 \r
1732 \r
1733  aodpid->SetITSsignal(track->GetITSsignal());\r
1734  aodpid->SetTPCsignal(track->GetTPCsignal());\r
1735  aodpid->SetTPCsignalN(track->GetTPCsignalN());\r
1736 \r
1737  //n TRD planes = 6\r
1738  Int_t nslices = track->GetNumberOfTRDslices()*6;\r
1739  Double_t *trdslices = new Double_t[nslices];\r
1740  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {\r
1741    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);\r
1742  }\r
1743  \r
1744 //TRD momentum\r
1745  for(Int_t iPl=0;iPl<6;iPl++){\r
1746    Double_t trdmom=track->GetTRDmomentum(iPl);\r
1747    aodpid->SetTRDmomentum(iPl,trdmom);\r
1748  }\r
1749 \r
1750  aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);\r
1751 \r
1752  //TOF PID  \r
1753  Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);\r
1754  aodpid->SetIntegratedTimes(times);\r
1755 \r
1756   Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());\r
1757   aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);\r
1758   \r
1759   Double_t tofRes[5];\r
1760   for (Int_t iMass=0; iMass<5; iMass++){\r
1761     tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));\r
1762   }\r
1763   aodpid->SetTOFpidResolution(tofRes);\r
1764 \r
1765   aodpid->SetHMPIDsignal(track->GetHMPIDsignal());\r
1766 \r
1767  //Extrapolate track to EMCAL surface for AOD-level track-cluster matching\r
1768  Double_t emcpos[3] = {0.,0.,0.};\r
1769  Double_t emcmom[3] = {0.,0.,0.};\r
1770  aodpid->SetEMCALPosition(emcpos);\r
1771  aodpid->SetEMCALMomentum(emcmom);\r
1772 \r
1773  AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();\r
1774  if(!outerparam) return;\r
1775 \r
1776  //To be replaced by call to AliEMCALGeoUtils when the class becomes available\r
1777  Bool_t okpos = outerparam->GetXYZ(emcpos);\r
1778  Bool_t okmom = outerparam->GetPxPyPz(emcmom);\r
1779  if(!(okpos && okmom)) return;\r
1780 \r
1781  aodpid->SetEMCALPosition(emcpos);\r
1782  aodpid->SetEMCALMomentum(emcmom);\r
1783 \r
1784 }\r
1785 \r
1786 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)\r
1787 {\r
1788     // Calculate chi2 per ndf for track\r
1789     Int_t  nClustersTPC = track->GetTPCNcls();\r
1790 \r
1791     if ( nClustersTPC > 5) {\r
1792        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));\r
1793     } else {\r
1794        return (-1.);\r
1795     }\r
1796  }\r
1797 \r
1798 \r
1799 //______________________________________________________________________________\r
1800 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)\r
1801 {\r
1802 // Terminate analysis\r
1803 //\r
1804     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");\r
1805 }\r
1806 \r
1807 //______________________________________________________________________________\r
1808 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){\r
1809 // Print MC info\r
1810   if(!pStack)return;\r
1811   label = TMath::Abs(label);\r
1812   TParticle *part = pStack->Particle(label);\r
1813   Printf("########################");\r
1814   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());\r
1815   part->Print();\r
1816   TParticle* mother = part;\r
1817   Int_t imo = part->GetFirstMother();\r
1818   Int_t nprim = pStack->GetNprimary();\r
1819   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {\r
1820   while((imo >= nprim)) {\r
1821     mother =  pStack->Particle(imo);\r
1822     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());\r
1823     mother->Print();\r
1824     imo =  mother->GetFirstMother();\r
1825   }\r
1826   Printf("########################");\r
1827 }\r