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