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