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->RelateToVertexTPC(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     Float_t ndf = track->GetTPCNcls()+1 - 5 ;\r
1002     if(ndf>0){\r
1003       aodTrack->SetChi2perNDF(track->GetTPCchi2()/ndf);\r
1004     }\r
1005     else{\r
1006       aodTrack->SetChi2perNDF(-1);\r
1007     }\r
1008     aodTrack->SetFlags(track->GetStatus());\r
1009     aodTrack->SetTPCPointsF(track->GetTPCNclsF());\r
1010 \r
1011     delete track;\r
1012   } // end of loop on tracks\r
1013   \r
1014 }\r
1015 \r
1016 //______________________________________________________________________________\r
1017 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)\r
1018 {\r
1019   // Tracks (primary and orphan)\r
1020 \r
1021   AliCodeTimerAuto("",0);\r
1022   \r
1023   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));\r
1024   \r
1025   const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
1026   Double_t p[3] = { 0. };\r
1027   Double_t pos[3] = { 0. };\r
1028   Double_t covTr[21] = { 0. };\r
1029   Double_t pid[10] = { 0. };\r
1030   AliAODTrack* aodTrack(0x0);\r
1031   AliAODPid* detpid(0x0);\r
1032   \r
1033   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) \r
1034   {\r
1035     if (fUsedTrack[nTrack]) continue;\r
1036     \r
1037     AliESDtrack *esdTrack = esd.GetTrack(nTrack);\r
1038     UInt_t selectInfo = 0;\r
1039     //\r
1040     // Track selection\r
1041     if (fTrackFilter) {\r
1042             selectInfo = fTrackFilter->IsSelected(esdTrack);\r
1043             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;\r
1044     }\r
1045     \r
1046     \r
1047     esdTrack->GetPxPyPz(p);\r
1048     esdTrack->GetXYZ(pos);\r
1049     esdTrack->GetCovarianceXYZPxPyPz(covTr);\r
1050     esdTrack->GetESDpid(pid);\r
1051     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());\r
1052     fPrimaryVertex->AddDaughter(aodTrack =\r
1053                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),\r
1054                                                             esdTrack->GetLabel(),\r
1055                                                             p,\r
1056                                                             kTRUE,\r
1057                                                             pos,\r
1058                                                             kFALSE,\r
1059                                                             covTr, \r
1060                                                             (Short_t)esdTrack->GetSign(),\r
1061                                                             esdTrack->GetITSClusterMap(), \r
1062                                                             pid,\r
1063                                                             fPrimaryVertex,\r
1064                                                             kTRUE, // check if this is right\r
1065                                                             vtx->UsesTrack(esdTrack->GetID()),\r
1066                                                             AliAODTrack::kPrimary, \r
1067                                                             selectInfo)\r
1068                          );\r
1069     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());\r
1070     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());\r
1071     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));\r
1072     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());\r
1073 \r
1074     fAODTrackRefs->AddAt(aodTrack, nTrack);\r
1075     \r
1076     \r
1077     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1078     aodTrack->SetFlags(esdTrack->GetStatus());\r
1079     aodTrack->ConvertAliPIDtoAODPID();\r
1080     SetAODPID(esdTrack,aodTrack,detpid);\r
1081   } // end of loop on tracks\r
1082 }\r
1083 \r
1084 //______________________________________________________________________________\r
1085 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)\r
1086 {\r
1087   AliCodeTimerAuto("",0);\r
1088   Int_t jPmdClusters=0;\r
1089   // Access to the AOD container of PMD clusters\r
1090   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());\r
1091   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {\r
1092     // file pmd clusters, to be revised!\r
1093     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);\r
1094     Int_t nLabel = 0;\r
1095     Int_t *label = 0x0;\r
1096     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};\r
1097     Double_t pidPmd[13] = { 0.}; // to be revised!\r
1098     // type not set!\r
1099     // assoc cluster not set\r
1100     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);\r
1101   }\r
1102 }\r
1103 \r
1104 //______________________________________________________________________________\r
1105 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)\r
1106 {\r
1107   AliCodeTimerAuto("",0);\r
1108   \r
1109   // Access to the AOD container of clusters\r
1110   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());\r
1111   Int_t jClusters(0);\r
1112   \r
1113   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {\r
1114     \r
1115     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);\r
1116     \r
1117     Int_t  id        = cluster->GetID();\r
1118     Int_t  nLabel    = cluster->GetNLabels();\r
1119     Int_t *labels    = cluster->GetLabels();\r
1120     if(labels){ \r
1121                   for(int i = 0;i < nLabel;++i){\r
1122                           if(fMChandler)fMChandler->SelectParticle(labels[i]);\r
1123                   }\r
1124           }             \r
1125     \r
1126     Float_t energy = cluster->E();\r
1127     Float_t posF[3] = { 0.};\r
1128     cluster->GetPosition(posF);\r
1129     \r
1130     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,\r
1131                                                                                       nLabel,\r
1132                                                                                       labels,\r
1133                                                                                       energy,\r
1134                                                                                       posF,\r
1135                                                                                       NULL,\r
1136                                                                                       cluster->GetType(),0);\r
1137     \r
1138     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),\r
1139                                 cluster->GetDispersion(),\r
1140                                 cluster->GetM20(), cluster->GetM02(),\r
1141                                 cluster->GetEmcCpvDistance(),  \r
1142                                 cluster->GetNExMax(),cluster->GetTOF()) ;\r
1143     \r
1144     caloCluster->SetPIDFromESD(cluster->GetPID());\r
1145     caloCluster->SetNCells(cluster->GetNCells());\r
1146     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());\r
1147     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());\r
1148     \r
1149     TArrayI* matchedT =         cluster->GetTracksMatched();\r
1150     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {        \r
1151       for (Int_t im = 0; im < matchedT->GetSize(); im++) {\r
1152         Int_t iESDtrack = matchedT->At(im);;\r
1153         if (fAODTrackRefs->At(iESDtrack) != 0) {\r
1154           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));\r
1155         }\r
1156       }\r
1157     }\r
1158     \r
1159   } \r
1160   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      \r
1161 }\r
1162 \r
1163 //______________________________________________________________________________\r
1164 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)\r
1165 {\r
1166   AliCodeTimerAuto("",0);\r
1167   // fill EMCAL cell info\r
1168   if (esd.GetEMCALCells()) { // protection against missing ESD information\r
1169     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());\r
1170     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;\r
1171     \r
1172     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());\r
1173     aodEMcells.CreateContainer(nEMcell);\r
1174     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);\r
1175     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      \r
1176       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell));\r
1177     }\r
1178     aodEMcells.Sort();\r
1179   }\r
1180 }\r
1181 \r
1182 //______________________________________________________________________________\r
1183 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)\r
1184 {\r
1185   AliCodeTimerAuto("",0);\r
1186   // fill PHOS cell info\r
1187   if (esd.GetPHOSCells()) { // protection against missing ESD information\r
1188     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());\r
1189     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;\r
1190     \r
1191     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());\r
1192     aodPHcells.CreateContainer(nPHcell);\r
1193     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);\r
1194     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      \r
1195       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));\r
1196     }\r
1197     aodPHcells.Sort();\r
1198   }\r
1199 }\r
1200 \r
1201 //______________________________________________________________________________\r
1202 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)\r
1203 {\r
1204   // tracklets    \r
1205   AliCodeTimerAuto("",0);\r
1206 \r
1207   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());\r
1208   const AliMultiplicity *mult = esd.GetMultiplicity();\r
1209   if (mult) {\r
1210     if (mult->GetNumberOfTracklets()>0) {\r
1211       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());\r
1212       \r
1213       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {\r
1214         if(fMChandler){\r
1215           fMChandler->SelectParticle(mult->GetLabel(n, 0));\r
1216           fMChandler->SelectParticle(mult->GetLabel(n, 1));\r
1217         }\r
1218         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));\r
1219       }\r
1220     }\r
1221   } else {\r
1222     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");\r
1223   }\r
1224 }\r
1225 \r
1226 //______________________________________________________________________________\r
1227 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)\r
1228 {\r
1229   AliCodeTimerAuto("",0);\r
1230   \r
1231   // Kinks: it is a big mess the access to the information in the kinks\r
1232   // The loop is on the tracks in order to find the mother and daugther of each kink\r
1233   \r
1234   Double_t covTr[21]={0.};\r
1235   Double_t pid[10]={0.};\r
1236   AliAODPid* detpid(0x0);\r
1237   \r
1238   fNumberOfKinks = esd.GetNumberOfKinks();\r
1239 \r
1240   const AliESDVertex* vtx = esd.GetPrimaryVertex();\r
1241   \r
1242   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack) \r
1243   {\r
1244     AliESDtrack * esdTrack = esd.GetTrack(iTrack);\r
1245     \r
1246     Int_t ikink = esdTrack->GetKinkIndex(0);\r
1247     \r
1248     if (ikink && fNumberOfKinks) {\r
1249             // Negative kink index: mother, positive: daughter\r
1250             \r
1251             // Search for the second track of the kink\r
1252             \r
1253             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {\r
1254         \r
1255         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);\r
1256         \r
1257         Int_t jkink = esdTrack1->GetKinkIndex(0);\r
1258         \r
1259         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {\r
1260           \r
1261           // The two tracks are from the same kink\r
1262           \r
1263           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks\r
1264           \r
1265           Int_t imother = -1;\r
1266           Int_t idaughter = -1;\r
1267           \r
1268           if (ikink<0 && jkink>0) {\r
1269             \r
1270             imother = iTrack;\r
1271             idaughter = jTrack;\r
1272           }\r
1273           else if (ikink>0 && jkink<0) {\r
1274             \r
1275             imother = jTrack;\r
1276             idaughter = iTrack;\r
1277           }\r
1278           else {\r
1279             //                  cerr << "Error: Wrong combination of kink indexes: "\r
1280             //                       << ikink << " " << jkink << endl;\r
1281             continue;\r
1282           }\r
1283           \r
1284           // Add the mother track if it passed primary track selection cuts\r
1285           \r
1286           AliAODTrack * mother = NULL;\r
1287           \r
1288           UInt_t selectInfo = 0;\r
1289           if (fTrackFilter) {\r
1290             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));\r
1291             if (!selectInfo) continue;\r
1292           }\r
1293           \r
1294           if (!fUsedTrack[imother]) {\r
1295             \r
1296             fUsedTrack[imother] = kTRUE;\r
1297             \r
1298             AliESDtrack *esdTrackM = esd.GetTrack(imother);\r
1299             Double_t p[3] = { 0. };\r
1300             Double_t pos[3] = { 0. };\r
1301             esdTrackM->GetPxPyPz(p);\r
1302             esdTrackM->GetXYZ(pos);\r
1303             esdTrackM->GetCovarianceXYZPxPyPz(covTr);\r
1304             esdTrackM->GetESDpid(pid);\r
1305             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());\r
1306             mother = \r
1307             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),\r
1308                                                esdTrackM->GetLabel(),\r
1309                                                p,\r
1310                                                kTRUE,\r
1311                                                pos,\r
1312                                                kFALSE,\r
1313                                                covTr, \r
1314                                                (Short_t)esdTrackM->GetSign(),\r
1315                                                esdTrackM->GetITSClusterMap(), \r
1316                                                pid,\r
1317                                                fPrimaryVertex,\r
1318                                                kTRUE, // check if this is right\r
1319                                                vtx->UsesTrack(esdTrack->GetID()),\r
1320                                                AliAODTrack::kPrimary,\r
1321                                                selectInfo);\r
1322             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());\r
1323             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());\r
1324             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));\r
1325             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());\r
1326 \r
1327             fAODTrackRefs->AddAt(mother, imother);\r
1328             \r
1329             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1330             mother->SetFlags(esdTrackM->GetStatus());\r
1331             mother->ConvertAliPIDtoAODPID();\r
1332             fPrimaryVertex->AddDaughter(mother);\r
1333             mother->ConvertAliPIDtoAODPID();\r
1334             SetAODPID(esdTrackM,mother,detpid);\r
1335           }\r
1336           else {\r
1337             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1\r
1338             //                       << " track " << imother << " has already been used!" << endl;\r
1339           }\r
1340           \r
1341           // Add the kink vertex\r
1342           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);\r
1343           \r
1344           AliAODVertex * vkink = \r
1345           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),\r
1346                                                   NULL,\r
1347                                                   0.,\r
1348                                                   mother,\r
1349                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!\r
1350                                                   AliAODVertex::kKink);\r
1351           // Add the daughter track\r
1352           \r
1353           AliAODTrack * daughter = NULL;\r
1354           \r
1355           if (!fUsedTrack[idaughter]) {\r
1356             \r
1357             fUsedTrack[idaughter] = kTRUE;\r
1358             \r
1359             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);\r
1360             Double_t p[3] = { 0. };\r
1361             Double_t pos[3] = { 0. };\r
1362 \r
1363             esdTrackD->GetPxPyPz(p);\r
1364             esdTrackD->GetXYZ(pos);\r
1365             esdTrackD->GetCovarianceXYZPxPyPz(covTr);\r
1366             esdTrackD->GetESDpid(pid);\r
1367             selectInfo = 0;\r
1368             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);\r
1369             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());\r
1370             daughter = \r
1371             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),\r
1372                                                esdTrackD->GetLabel(),\r
1373                                                p,\r
1374                                                kTRUE,\r
1375                                                pos,\r
1376                                                kFALSE,\r
1377                                                covTr, \r
1378                                                (Short_t)esdTrackD->GetSign(),\r
1379                                                esdTrackD->GetITSClusterMap(), \r
1380                                                pid,\r
1381                                                vkink,\r
1382                                                kTRUE, // check if this is right\r
1383                                                vtx->UsesTrack(esdTrack->GetID()),\r
1384                                                AliAODTrack::kSecondary,\r
1385                                                selectInfo);\r
1386             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());\r
1387             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());\r
1388             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());\r
1389             fAODTrackRefs->AddAt(daughter, idaughter);\r
1390             \r
1391             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;\r
1392             daughter->SetFlags(esdTrackD->GetStatus());\r
1393             daughter->ConvertAliPIDtoAODPID();\r
1394             vkink->AddDaughter(daughter);\r
1395             daughter->ConvertAliPIDtoAODPID();\r
1396             SetAODPID(esdTrackD,daughter,detpid);\r
1397           }\r
1398           else {\r
1399             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1\r
1400             //                       << " track " << idaughter << " has already been used!" << endl;\r
1401           }\r
1402         }\r
1403             }\r
1404     }      \r
1405   }\r
1406 }\r
1407 \r
1408 //______________________________________________________________________________\r
1409 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)\r
1410 {\r
1411   AliCodeTimerAuto("",0);\r
1412   \r
1413   // Access to the AOD container of vertices\r
1414   fNumberOfVertices = 0;\r
1415   \r
1416   Double_t pos[3] = { 0. };\r
1417   Double_t covVtx[6] = { 0. };\r
1418 \r
1419   // Add primary vertex. The primary tracks will be defined\r
1420   // after the loops on the composite objects (V0, cascades, kinks)\r
1421   const AliESDVertex *vtx = esd.GetPrimaryVertex();\r
1422   \r
1423   vtx->GetXYZ(pos); // position\r
1424   vtx->GetCovMatrix(covVtx); //covariance matrix\r
1425   \r
1426   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])\r
1427   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);\r
1428   fPrimaryVertex->SetName(vtx->GetName());\r
1429   fPrimaryVertex->SetTitle(vtx->GetTitle());\r
1430   \r
1431   TString vtitle = vtx->GetTitle();\r
1432   if (!vtitle.Contains("VertexerTracks")) \r
1433     fPrimaryVertex->SetNContributors(vtx->GetNContributors());\r
1434   \r
1435   if (fDebug > 0) fPrimaryVertex->Print();  \r
1436   \r
1437   // Add SPD "main" vertex \r
1438   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();\r
1439   vtxS->GetXYZ(pos); // position\r
1440   vtxS->GetCovMatrix(covVtx); //covariance matrix\r
1441   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])\r
1442   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);\r
1443   mVSPD->SetName(vtxS->GetName());\r
1444   mVSPD->SetTitle(vtxS->GetTitle());\r
1445   mVSPD->SetNContributors(vtxS->GetNContributors()); \r
1446   \r
1447   // Add SPD pileup vertices\r
1448   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)\r
1449   {\r
1450     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);\r
1451     vtxP->GetXYZ(pos); // position\r
1452     vtxP->GetCovMatrix(covVtx); //covariance matrix\r
1453     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])\r
1454     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);\r
1455     pVSPD->SetName(vtxP->GetName());\r
1456     pVSPD->SetTitle(vtxP->GetTitle());\r
1457     pVSPD->SetNContributors(vtxP->GetNContributors()); \r
1458   }\r
1459   \r
1460   // Add TRK pileup vertices\r
1461   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)\r
1462   {\r
1463     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);\r
1464     vtxP->GetXYZ(pos); // position\r
1465     vtxP->GetCovMatrix(covVtx); //covariance matrix\r
1466     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])\r
1467     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);\r
1468     pVTRK->SetName(vtxP->GetName());\r
1469     pVTRK->SetTitle(vtxP->GetTitle());\r
1470     pVTRK->SetNContributors(vtxP->GetNContributors());\r
1471   }\r
1472 }\r
1473 \r
1474 //______________________________________________________________________________\r
1475 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)\r
1476 {\r
1477   // Convert VZERO data\r
1478   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();\r
1479   *vzeroData = *(esd.GetVZEROData());\r
1480 }\r
1481 \r
1482 //______________________________________________________________________________\r
1483 void AliAnalysisTaskESDfilter::ConvertESDtoAOD() \r
1484 {\r
1485   // ESD Filter analysis task executed for each event\r
1486   \r
1487   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());\r
1488   \r
1489   if(!esd)return;\r
1490 \r
1491   AliCodeTimerAuto("",0);\r
1492   \r
1493   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );\r
1494   \r
1495   fNumberOfTracks = 0;\r
1496   fNumberOfPositiveTracks = 0;\r
1497   fNumberOfV0s = 0;\r
1498   fNumberOfVertices = 0;\r
1499   fNumberOfCascades = 0;\r
1500   fNumberOfKinks = 0;\r
1501     \r
1502   AliAODHeader* header = ConvertHeader(*esd);\r
1503 \r
1504   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);\r
1505   \r
1506   // Fetch Stack for debuggging if available \r
1507   fMChandler=0x0;\r
1508   if(MCEvent())\r
1509   {\r
1510     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); \r
1511   }\r
1512   \r
1513   // loop over events and fill them\r
1514   // Multiplicity information needed by the header (to be revised!)\r
1515   Int_t nTracks    = esd->GetNumberOfTracks();\r
1516   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);\r
1517 \r
1518   // Update the header\r
1519 \r
1520   Int_t nV0s      = esd->GetNumberOfV0s();\r
1521   Int_t nCascades = esd->GetNumberOfCascades();\r
1522   Int_t nKinks    = esd->GetNumberOfKinks();\r
1523   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;\r
1524   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex\r
1525   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();\r
1526   nVertices+=nPileSPDVertices;\r
1527   nVertices+=nPileTrkVertices;\r
1528   Int_t nJets     = 0;\r
1529   Int_t nCaloClus = esd->GetNumberOfCaloClusters();\r
1530   Int_t nFmdClus  = 0;\r
1531   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();\r
1532     \r
1533   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));\r
1534        \r
1535   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus);\r
1536 \r
1537   if (nV0s > 0) \r
1538   {\r
1539     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0\r
1540     fAODV0VtxRefs = new TRefArray(nV0s);\r
1541     // RefArray to store the mapping between esd V0 number and newly created AOD-V0\r
1542     fAODV0Refs = new TRefArray(nV0s); \r
1543     // Array to take into account the V0s already added to the AOD (V0 within cascades)\r
1544     fUsedV0 = new Bool_t[nV0s];\r
1545     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;\r
1546   }\r
1547   \r
1548   if (nTracks>0) \r
1549   {\r
1550     // RefArray to store the mapping between esd track number and newly created AOD-Track\r
1551     \r
1552     fAODTrackRefs = new TRefArray(nTracks);\r
1553 \r
1554     // Array to take into account the tracks already added to the AOD    \r
1555     fUsedTrack = new Bool_t[nTracks];\r
1556     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;\r
1557   }\r
1558   \r
1559   // Array to take into account the kinks already added to the AOD\r
1560   if (nKinks>0) \r
1561   {\r
1562     fUsedKink = new Bool_t[nKinks];\r
1563     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;\r
1564   }\r
1565     \r
1566   ConvertPrimaryVertices(*esd);\r
1567 \r
1568   //setting best TOF PID\r
1569   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
1570   if (esdH)\r
1571       fESDpid = esdH->GetESDpid();\r
1572 \r
1573   if (fIsPidOwner && fESDpid){\r
1574     delete fESDpid;\r
1575     fESDpid = 0;\r
1576   }\r
1577   if(!fESDpid)\r
1578   { //in case of no Tender attached \r
1579     fESDpid = new AliESDpid;\r
1580     fIsPidOwner = kTRUE;\r
1581   }\r
1582   \r
1583   if(!esd->GetTOFHeader())\r
1584   { //protection in case the pass2 LHC10b,c,d have been processed without tender. \r
1585     Float_t t0spread[10];\r
1586     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!! \r
1587     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
1588     fESDpid->GetTOFResponse().SetT0resolution(t0spread);\r
1589     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);\r
1590     \r
1591     fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);    \r
1592   }\r
1593   \r
1594   if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender. \r
1595   \r
1596   if ( fAreCascadesEnabled ) ConvertCascades(*esd);\r
1597 \r
1598   if ( fAreV0sEnabled ) ConvertV0s(*esd);\r
1599   \r
1600   if ( fAreKinksEnabled ) ConvertKinks(*esd);\r
1601   \r
1602   if ( fAreTracksEnabled ) ConvertTracks(*esd);\r
1603   \r
1604   // Update number of AOD tracks in header at the end of track loop (M.G.)\r
1605   header->SetRefMultiplicity(fNumberOfTracks);\r
1606   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);\r
1607   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);\r
1608 \r
1609   if ( fTPCOnlyFilterMask ) ConvertTPCOnlyTracks(*esd);\r
1610   \r
1611   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);\r
1612   \r
1613   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);\r
1614   \r
1615   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);\r
1616   \r
1617   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);\r
1618   \r
1619   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);\r
1620   \r
1621   delete fAODTrackRefs; fAODTrackRefs=0x0;\r
1622   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;\r
1623   delete fAODV0Refs; fAODV0Refs=0x0;\r
1624   \r
1625   delete[] fUsedTrack; fUsedTrack=0x0;\r
1626   delete[] fUsedV0; fUsedV0=0x0;\r
1627   delete[] fUsedKink; fUsedKink=0x0;\r
1628 \r
1629   if ( fIsPidOwner){\r
1630     delete fESDpid;\r
1631     fESDpid = 0x0;\r
1632   }\r
1633 \r
1634 \r
1635 }\r
1636 \r
1637 \r
1638 //______________________________________________________________________________\r
1639 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)\r
1640 {\r
1641   //\r
1642   // Setter for the raw PID detector signals\r
1643   //\r
1644 \r
1645   // Save PID object for candidate electrons\r
1646     Bool_t pidSave = kFALSE;\r
1647     if (fTrackFilter) {\r
1648         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");\r
1649         if (selectInfo)  pidSave = kTRUE;\r
1650     }\r
1651 \r
1652 \r
1653     // Tracks passing pt cut \r
1654     if(esdtrack->Pt()>fHighPthreshold) {\r
1655         pidSave = kTRUE;\r
1656     } else {\r
1657         if(fPtshape){\r
1658             if(esdtrack->Pt()> fPtshape->GetXmin()){\r
1659                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);\r
1660                 if(gRandom->Rndm(0)<1./y){\r
1661                     pidSave = kTRUE;\r
1662                 }//end rndm\r
1663             }//end if p < pmin\r
1664         }//end if p function\r
1665     }// end else\r
1666 \r
1667     if (pidSave) {\r
1668       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track\r
1669         detpid = new AliAODPid();\r
1670         SetDetectorRawSignals(detpid,esdtrack);\r
1671         aodtrack->SetDetPID(detpid);\r
1672       }\r
1673     }\r
1674 }\r
1675 \r
1676 //______________________________________________________________________________\r
1677 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)\r
1678 {\r
1679 //\r
1680 //assignment of the detector signals (AliXXXesdPID inspired)\r
1681 //\r
1682  if(!track){\r
1683  AliInfo("no ESD track found. .....exiting");\r
1684  return;\r
1685  }\r
1686  // TPC momentum\r
1687  const AliExternalTrackParam *in=track->GetInnerParam();\r
1688  if (in) {\r
1689    aodpid->SetTPCmomentum(in->GetP());\r
1690  }else{\r
1691    aodpid->SetTPCmomentum(-1.);\r
1692  }\r
1693 \r
1694 \r
1695  aodpid->SetITSsignal(track->GetITSsignal());\r
1696  aodpid->SetTPCsignal(track->GetTPCsignal());\r
1697  aodpid->SetTPCsignalN(track->GetTPCsignalN());\r
1698 \r
1699  //n TRD planes = 6\r
1700  Int_t nslices = track->GetNumberOfTRDslices()*6;\r
1701  Double_t *trdslices = new Double_t[nslices];\r
1702  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {\r
1703    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);\r
1704  }\r
1705  \r
1706 //TRD momentum\r
1707  for(Int_t iPl=0;iPl<6;iPl++){\r
1708    Double_t trdmom=track->GetTRDmomentum(iPl);\r
1709    aodpid->SetTRDmomentum(iPl,trdmom);\r
1710  }\r
1711 \r
1712  aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices);\r
1713 \r
1714  //TOF PID  \r
1715  Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);\r
1716  aodpid->SetIntegratedTimes(times);\r
1717 \r
1718   Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());\r
1719   aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);\r
1720   \r
1721   Double_t tofRes[5];\r
1722   for (Int_t iMass=0; iMass<5; iMass++){\r
1723     tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));\r
1724   }\r
1725   aodpid->SetTOFpidResolution(tofRes);\r
1726 \r
1727   aodpid->SetHMPIDsignal(track->GetHMPIDsignal());\r
1728 \r
1729  //Extrapolate track to EMCAL surface for AOD-level track-cluster matching\r
1730  Double_t emcpos[3] = {0.,0.,0.};\r
1731  Double_t emcmom[3] = {0.,0.,0.};\r
1732  aodpid->SetEMCALPosition(emcpos);\r
1733  aodpid->SetEMCALMomentum(emcmom);\r
1734 \r
1735  AliExternalTrackParam *outerparam = (AliExternalTrackParam*)track->GetOuterParam();\r
1736  if(!outerparam) return;\r
1737 \r
1738  //To be replaced by call to AliEMCALGeoUtils when the class becomes available\r
1739  Bool_t okpos = outerparam->GetXYZ(emcpos);\r
1740  Bool_t okmom = outerparam->GetPxPyPz(emcmom);\r
1741  if(!(okpos && okmom)) return;\r
1742 \r
1743  aodpid->SetEMCALPosition(emcpos);\r
1744  aodpid->SetEMCALMomentum(emcmom);\r
1745 \r
1746 }\r
1747 \r
1748 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)\r
1749 {\r
1750     // Calculate chi2 per ndf for track\r
1751     Int_t  nClustersTPC = track->GetTPCNcls();\r
1752 \r
1753     if ( nClustersTPC > 5) {\r
1754        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));\r
1755     } else {\r
1756        return (-1.);\r
1757     }\r
1758  }\r
1759 \r
1760 \r
1761 //______________________________________________________________________________\r
1762 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)\r
1763 {\r
1764 // Terminate analysis\r
1765 //\r
1766     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");\r
1767 }\r
1768 \r
1769 //______________________________________________________________________________\r
1770 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){\r
1771 // Print MC info\r
1772   if(!pStack)return;\r
1773   label = TMath::Abs(label);\r
1774   TParticle *part = pStack->Particle(label);\r
1775   Printf("########################");\r
1776   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());\r
1777   part->Print();\r
1778   TParticle* mother = part;\r
1779   Int_t imo = part->GetFirstMother();\r
1780   Int_t nprim = pStack->GetNprimary();\r
1781   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {\r
1782   while((imo >= nprim)) {\r
1783     mother =  pStack->Particle(imo);\r
1784     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());\r
1785     mother->Print();\r
1786     imo =  mother->GetFirstMother();\r
1787   }\r
1788   Printf("########################");\r
1789 }\r