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