ANALYSIS/AliAnalysisTaskESDfilter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 #include <TChain.h>
  19 #include <TTree.h>
  20 #include <TList.h>
  21 #include <TArrayI.h>
  22 #include <TParameter.h>
  23 #include <TRandom.h>
  24 #include <TParticle.h>
  25 #include <TFile.h>
  26
  27 #include "AliAnalysisTaskESDfilter.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliESDEvent.h"
  30 #include "AliESDRun.h"
  31 #include "AliStack.h"
  32 #include "AliAODEvent.h"
  33 #include "AliMCEvent.h"
  34 #include "AliMCEventHandler.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAODMCParticle.h"
  38 #include "AliAnalysisFilter.h"
  39 #include "AliESDMuonTrack.h"
  40 #include "AliESDVertex.h"
  41 #include "AliCentrality.h"
  42 #include "AliEventplane.h"
  43 #include "AliESDv0.h"
  44 #include "AliESDkink.h"
  45 #include "AliESDcascade.h"
  46 #include "AliESDPmdTrack.h"
  47 #include "AliESDCaloCluster.h"
  48 #include "AliESDCaloCells.h"
  49 #include "AliMultiplicity.h"
  50 #include "AliLog.h"
  51 #include "AliCodeTimer.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "AliESDpid.h"
  54 #include "AliAODHMPIDrings.h"
  55 #include "AliV0vertexer.h"
  56 #include "AliCascadeVertexer.h"
  57 #include "Riostream.h"
  58 #include "AliExternalTrackParam.h"
  59 #include "AliTrackerBase.h"
  60 #include "TVector3.h"
  61 #include "AliTPCdEdxInfo.h"
  62
  63 using std::cout;
  64 using std::endl;
  65 ClassImp(AliAnalysisTaskESDfilter)
  66
  67 ////////////////////////////////////////////////////////////////////////
  68
  69 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():
  70   AliAnalysisTaskSE(),
  71   fTrackFilter(0x0),
  72   fKinkFilter(0x0),
  73   fV0Filter(0x0),
  74   fCascadeFilter(0x0),
  75   fHighPthreshold(0),
  76   fPtshape(0x0),
  77   fEnableFillAOD(kTRUE),
  78   fUsedTrack(0x0),
  79   fUsedKink(0x0),
  80   fUsedV0(0x0),
  81   fAODTrackRefs(0x0),
  82   fAODV0VtxRefs(0x0),
  83   fAODV0Refs(0x0),
  84   fMChandler(0x0),
  85   fNumberOfTracks(0),
  86   fNumberOfPositiveTracks(0),
  87   fNumberOfV0s(0),
  88   fNumberOfVertices(0),
  89   fNumberOfCascades(0),
  90   fNumberOfKinks(0),
  91   fOldESDformat(kFALSE),
  92   fPrimaryVertex(0x0),
  93   fTPCConstrainedFilterMask(0),
  94   fHybridFilterMaskTPCCG(0),
  95   fWriteHybridTPCCOnly(kFALSE),
  96   fGlobalConstrainedFilterMask(0),
  97   fHybridFilterMaskGCG(0),
  98   fWriteHybridGCOnly(kFALSE),
  99   fIsVZEROEnabled(kTRUE),
 100   fIsTZEROEnabled(kTRUE),
 101   fIsZDCEnabled(kTRUE),
 102   fIsHMPIDEnabled(kTRUE),
 103   fIsV0CascadeRecoEnabled(kFALSE),
 104   fAreCascadesEnabled(kTRUE),
 105   fAreV0sEnabled(kTRUE),
 106   fAreKinksEnabled(kTRUE),
 107   fAreTracksEnabled(kTRUE),
 108   fArePmdClustersEnabled(kTRUE),
 109   fAreCaloClustersEnabled(kTRUE),
 110   fAreEMCALCellsEnabled(kTRUE),
 111   fArePHOSCellsEnabled(kTRUE),
 112   fAreEMCALTriggerEnabled(kTRUE),
 113   fArePHOSTriggerEnabled(kTRUE),
 114   fAreTrackletsEnabled(kTRUE),
 115   fESDpid(0x0),
 116   fIsPidOwner(kFALSE),
 117   fTPCaloneTrackCuts(0),
 118   fDoPropagateTrackToEMCal(kTRUE)
 119 {
 120   // Default constructor
 121     fV0Cuts[0] =  33.   ;   // max allowed chi2
 122     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 123     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 124     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 125     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 126     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 127     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 128
 129     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 130     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 131     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 132     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 133     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 134     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 135     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 136     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 137 }
 138
 139 //______________________________________________________________________________
 140 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):
 141     AliAnalysisTaskSE(name),
 142     fTrackFilter(0x0),
 143     fKinkFilter(0x0),
 144     fV0Filter(0x0),
 145     fCascadeFilter(0x0),
 146     fHighPthreshold(0),
 147     fPtshape(0x0),
 148     fEnableFillAOD(kTRUE),
 149     fUsedTrack(0x0),
 150     fUsedKink(0x0),
 151     fUsedV0(0x0),
 152     fAODTrackRefs(0x0),
 153     fAODV0VtxRefs(0x0),
 154     fAODV0Refs(0x0),
 155     fMChandler(0x0),
 156     fNumberOfTracks(0),
 157     fNumberOfPositiveTracks(0),
 158     fNumberOfV0s(0),
 159     fNumberOfVertices(0),
 160     fNumberOfCascades(0),
 161     fNumberOfKinks(0),
 162     fOldESDformat(kFALSE),
 163     fPrimaryVertex(0x0),
 164   fTPCConstrainedFilterMask(0),
 165   fHybridFilterMaskTPCCG(0),
 166   fWriteHybridTPCCOnly(kFALSE),
 167   fGlobalConstrainedFilterMask(0),
 168   fHybridFilterMaskGCG(0),
 169   fWriteHybridGCOnly(kFALSE),
 170     fIsVZEROEnabled(kTRUE),
 171     fIsTZEROEnabled(kTRUE),
 172     fIsZDCEnabled(kTRUE),
 173     fIsHMPIDEnabled(kTRUE),
 174     fIsV0CascadeRecoEnabled(kFALSE),
 175     fAreCascadesEnabled(kTRUE),
 176     fAreV0sEnabled(kTRUE),
 177     fAreKinksEnabled(kTRUE),
 178     fAreTracksEnabled(kTRUE),
 179     fArePmdClustersEnabled(kTRUE),
 180     fAreCaloClustersEnabled(kTRUE),
 181     fAreEMCALCellsEnabled(kTRUE),
 182     fArePHOSCellsEnabled(kTRUE),
 183                 fAreEMCALTriggerEnabled(kTRUE),
 184                 fArePHOSTriggerEnabled(kTRUE),
 185                 fAreTrackletsEnabled(kTRUE),
 186     fESDpid(0x0),
 187     fIsPidOwner(kFALSE),
 188     fTPCaloneTrackCuts(0),
 189   fDoPropagateTrackToEMCal(kTRUE)
 190 {
 191   // Constructor
 192
 193     fV0Cuts[0] =  33.   ;   // max allowed chi2
 194     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 195     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 196     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 197     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 198     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 199     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 200
 201     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 202     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 203     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 204     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 205     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 206     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 207     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 208     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 209
 210
 211
 212 }
 213 AliAnalysisTaskESDfilter::~AliAnalysisTaskESDfilter(){
 214     if(fIsPidOwner) delete fESDpid;
 215 }
 216 //______________________________________________________________________________
 217 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()
 218 {
 219   //
 220   // Create Output Objects conenct filter to outputtree
 221   //
 222   if(OutputTree())
 223   {
 224     OutputTree()->GetUserInfo()->Add(fTrackFilter);
 225   }
 226   else
 227   {
 228     AliError("No OutputTree() for adding the track filter");
 229   }
 230   fTPCaloneTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 231 }
 232
 233 //______________________________________________________________________________
 234 void AliAnalysisTaskESDfilter::Init()
 235 {
 236   // Initialization
 237   if (fDebug > 1) AliInfo("Init() \n");
 238   // Call configuration file
 239 }
 240
 241 //______________________________________________________________________________
 242 Bool_t AliAnalysisTaskESDfilter::Notify()
 243 {
 244 // Notify method.
 245    AddMetadataToUserInfo();
 246    return kTRUE;
 247 }
 248
 249 //______________________________________________________________________________
 250 Bool_t AliAnalysisTaskESDfilter::AddMetadataToUserInfo()
 251 {
 252 // Copy metadata to AOD user info.
 253   static Bool_t copyFirst = kFALSE;
 254   if (!copyFirst) {
 255     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 256     if (!mgr) {
 257        AliError("AliAnalysisTaskESDfilter::AddMetadataToUserInfo() : No analysis manager !");
 258        return kFALSE;
 259     }
 260     TTree *esdTree = mgr->GetTree()->GetTree();
 261     if (!esdTree) return kFALSE;
 262     TNamed *alirootVersion = (TNamed*)esdTree->GetUserInfo()->FindObject("alirootVersion");
 263     if (!alirootVersion) return kFALSE;
 264     AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
 265          if (!aodHandler) return kFALSE;
 266     TTree *aodTree = aodHandler->GetTree();
 267     if (!aodTree) return kFALSE;
 268     aodTree->GetUserInfo()->Add(new TNamed(*alirootVersion));
 269     copyFirst = kTRUE;
 270   }
 271   return kTRUE;
 272 }
 273
 274 //______________________________________________________________________________
 275 void AliAnalysisTaskESDfilter::PrintTask(Option_t *option, Int_t indent) const
 276 {
 277 // Print selection task information
 278   AliInfo("");
 279
 280   AliAnalysisTaskSE::PrintTask(option,indent);
 281
 282   TString spaces(' ',indent+3);
 283
 284         cout << spaces.Data() << Form("Cascades       are %s",fAreCascadesEnabled ? "ENABLED":"DISABLED") << endl;
 285         cout << spaces.Data() << Form("V0s            are %s",fAreV0sEnabled ? "ENABLED":"DISABLED") << endl;
 286         cout << spaces.Data() << Form("Kinks          are %s",fAreKinksEnabled ? "ENABLED":"DISABLED") << endl;
 287         cout << spaces.Data() << Form("Tracks         are %s",fAreTracksEnabled ? "ENABLED":"DISABLED") << endl;
 288         cout << spaces.Data() << Form("PmdClusters    are %s",fArePmdClustersEnabled ? "ENABLED":"DISABLED") << endl;
 289         cout << spaces.Data() << Form("CaloClusters   are %s",fAreCaloClustersEnabled ? "ENABLED":"DISABLED") << endl;
 290   cout << spaces.Data() << Form("EMCAL cells    are %s",fAreEMCALCellsEnabled ? "ENABLED":"DISABLED") << endl;
 291         cout << spaces.Data() << Form("EMCAL triggers are %s",fAreEMCALTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 292         cout << spaces.Data() << Form("PHOS triggers  are %s",fArePHOSTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 293         cout << spaces.Data() << Form("Tracklets      are %s",fAreTrackletsEnabled ? "ENABLED":"DISABLED") << endl;
 294         cout << spaces.Data() << Form("PropagateTrackToEMCal  is %s", fDoPropagateTrackToEMCal ? "ENABLED":"DISABLED") << endl;
 295 }
 296
 297 //______________________________________________________________________________
 298 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)
 299 {
 300 // Execute analysis for current event
 301 //
 302
 303   Long64_t ientry = Entry();
 304
 305   if (fDebug > 0) {
 306       printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
 307       if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
 308       if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
 309   }
 310   // Filters must explicitely enable AOD filling in their UserExec (AG)
 311   if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
 312   if(fEnableFillAOD) {
 313      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 314      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 315   }
 316   ConvertESDtoAOD();
 317 }
 318
 319 //______________________________________________________________________________
 320 TClonesArray& AliAnalysisTaskESDfilter::Cascades()
 321 {
 322   return *(AODEvent()->GetCascades());
 323 }
 324
 325 //______________________________________________________________________________
 326 TClonesArray& AliAnalysisTaskESDfilter::Tracks()
 327 {
 328   return *(AODEvent()->GetTracks());
 329 }
 330
 331 //______________________________________________________________________________
 332 TClonesArray& AliAnalysisTaskESDfilter::V0s()
 333 {
 334   return *(AODEvent()->GetV0s());
 335 }
 336
 337 //______________________________________________________________________________
 338 TClonesArray& AliAnalysisTaskESDfilter::Vertices()
 339 {
 340   return *(AODEvent()->GetVertices());
 341 }
 342
 343 //______________________________________________________________________________
 344 AliAODHeader* AliAnalysisTaskESDfilter::ConvertHeader(const AliESDEvent& esd)
 345 {
 346 // Convert header information
 347
 348   AliCodeTimerAuto("",0);
 349
 350   AliAODHeader* header = AODEvent()->GetHeader();
 351
 352   header->SetRunNumber(esd.GetRunNumber());
 353   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 354
 355   TTree* tree = fInputHandler->GetTree();
 356   if (tree) {
 357     TFile* file = tree->GetCurrentFile();
 358     if (file) header->SetESDFileName(file->GetName());
 359   }
 360
 361   if (fOldESDformat) {
 362     header->SetBunchCrossNumber(0);
 363     header->SetOrbitNumber(0);
 364     header->SetPeriodNumber(0);
 365     header->SetEventType(0);
 366     header->SetMuonMagFieldScale(-999.);
 367     header->SetCentrality(0);
 368     header->SetEventplane(0);
 369   } else {
 370     header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
 371     header->SetOrbitNumber(esd.GetOrbitNumber());
 372     header->SetPeriodNumber(esd.GetPeriodNumber());
 373     header->SetEventType(esd.GetEventType());
 374
 375     header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
 376     if(const_cast<AliESDEvent&>(esd).GetCentrality()){
 377       header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
 378     }
 379     else{
 380       header->SetCentrality(0);
 381     }
 382     if(const_cast<AliESDEvent&>(esd).GetEventplane()){
 383       header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
 384     }
 385     else{
 386       header->SetEventplane(0);
 387     }
 388   }
 389
 390   // Trigger
 391   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
 392   header->SetTriggerMask(esd.GetTriggerMask());
 393   header->SetTriggerCluster(esd.GetTriggerCluster());
 394   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());
 395   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());
 396   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());
 397
 398   header->SetMagneticField(esd.GetMagneticField());
 399   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
 400   header->SetZDCN1Energy(esd.GetZDCN1Energy());
 401   header->SetZDCP1Energy(esd.GetZDCP1Energy());
 402   header->SetZDCN2Energy(esd.GetZDCN2Energy());
 403   header->SetZDCP2Energy(esd.GetZDCP2Energy());
 404   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
 405
 406   // ITS Cluster Multiplicty
 407   const AliMultiplicity *mult = esd.GetMultiplicity();
 408   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
 409
 410   // TPC only Reference Multiplicty
 411   Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
 412   header->SetTPConlyRefMultiplicity(refMult);
 413
 414   //
 415   Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};
 416   Float_t diamcov[3];
 417   esd.GetDiamondCovXY(diamcov);
 418   header->SetDiamond(diamxy,diamcov);
 419   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
 420
 421   // VZERO channel equalization factors for event-plane reconstruction
 422   header->SetVZEROEqFactors(esd.GetVZEROEqFactors());
 423
 424   // T0 Resolution information
 425   const AliESDRun* esdRun = esd.GetESDRun();
 426   for (Int_t i=0;i<AliESDRun::kT0spreadSize;i++) header->SetT0spread(i,esdRun->GetT0spread(i));
 427
 428   return header;
 429 }
 430
 431 //______________________________________________________________________________
 432 void AliAnalysisTaskESDfilter::ConvertCascades(const AliESDEvent& esd)
 433 {
 434
 435   // Convert the cascades part of the ESD.
 436   // Return the number of cascades
 437
 438   AliCodeTimerAuto("",0);
 439
 440   // Create vertices starting from the most complex objects
 441   Double_t chi2 = 0.;
 442
 443   const AliESDVertex* vtx = esd.GetPrimaryVertex();
 444   Double_t pos[3] = { 0. };
 445   Double_t covVtx[6] = { 0. };
 446   Double_t momBach[3]={0.};
 447   Double_t covTr[21]={0.};
 448   Double_t pid[10]={0.};
 449   AliAODPid* detpid(0x0);
 450   AliAODVertex* vV0FromCascade(0x0);
 451   AliAODv0* aodV0(0x0);
 452   AliAODcascade* aodCascade(0x0);
 453   AliAODTrack* aodTrack(0x0);
 454   Double_t momPos[3]={0.};
 455   Double_t momNeg[3] = { 0. };
 456   Double_t momPosAtV0vtx[3]={0.};
 457   Double_t momNegAtV0vtx[3]={0.};
 458
 459   TClonesArray& verticesArray = Vertices();
 460   TClonesArray& tracksArray = Tracks();
 461   TClonesArray& cascadesArray = Cascades();
 462
 463   // Cascades (Modified by A.Maire - February 2009)
 464   for (Int_t nCascade = 0; nCascade < esd.GetNumberOfCascades(); ++nCascade) {
 465
 466     // 0- Preparation
 467     //
 468     AliESDcascade *esdCascade = esd.GetCascade(nCascade);
 469                 Int_t  idxPosFromV0Dghter  = esdCascade->GetPindex();
 470                 Int_t  idxNegFromV0Dghter  = esdCascade->GetNindex();
 471                 Int_t  idxBachFromCascade  = esdCascade->GetBindex();
 472
 473     AliESDtrack  *esdCascadePos  = esd.GetTrack( idxPosFromV0Dghter);
 474     AliESDtrack  *esdCascadeNeg  = esd.GetTrack( idxNegFromV0Dghter);
 475     AliESDtrack  *esdCascadeBach = esd.GetTrack( idxBachFromCascade);
 476
 477     // Identification of the V0 within the esdCascade (via both daughter track indices)
 478     AliESDv0 * currentV0   = 0x0;
 479     Int_t      idxV0FromCascade = -1;
 480
 481     for (Int_t iV0=0; iV0<esd.GetNumberOfV0s(); ++iV0) {
 482
 483       currentV0 = esd.GetV0(iV0);
 484       Int_t posCurrentV0 = currentV0->GetPindex();
 485       Int_t negCurrentV0 = currentV0->GetNindex();
 486
 487       if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {
 488         idxV0FromCascade = iV0;
 489         break;
 490       }
 491     }
 492
 493     if(idxV0FromCascade < 0){
 494       printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");
 495       continue;
 496     }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"
 497
 498     AliESDv0 *esdV0FromCascade   = esd.GetV0(idxV0FromCascade);
 499
 500     // 1 - Cascade selection
 501
 502     //  AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd.GetPrimaryVertex()));
 503     //  TList cascadeObjects;
 504     //  cascadeObjects.AddAt(esdV0FromCascade, 0);
 505     //  cascadeObjects.AddAt(esdCascadePos,    1);
 506     //  cascadeObjects.AddAt(esdCascadeNeg,    2);
 507     //  cascadeObjects.AddAt(esdCascade,       3);
 508     //  cascadeObjects.AddAt(esdCascadeBach,   4);
 509     //  cascadeObjects.AddAt(esdPrimVtx,       5);
 510     //
 511     //  UInt_t selectCascade = 0;
 512     //  if (fCascadeFilter) {
 513     //    // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects);
 514     //          // FIXME AliESDCascadeCuts to be implemented ...
 515     //
 516     //          // Here we may encounter a moot point at the V0 level
 517     //          // between the cascade selections and the V0 ones :
 518     //          // the V0 selected along with the cascade (secondary V0) may
 519     //          // usually be removed from the dedicated V0 selections (prim V0) ...
 520     //          // -> To be discussed !
 521     //
 522     //    // this is a little awkward but otherwise the
 523     //    // list wants to access the pointer (delete it)
 524     //    // again when going out of scope
 525     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 526     //    esdPrimVtx = 0;
 527     //    if (!selectCascade)
 528     //      continue;
 529     //  }
 530     //  else{
 531     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 532     //    esdPrimVtx = 0;
 533     //  }
 534
 535     // 2 - Add the cascade vertex
 536
 537     esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);
 538     esdCascade->GetPosCovXi(covVtx);
 539     chi2 = esdCascade->GetChi2Xi();
 540
 541     AliAODVertex *vCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex( pos,
 542                                                                      covVtx,
 543                                                                      chi2, // FIXME = Chi2/NDF will be needed
 544                                                                      fPrimaryVertex,
 545                                                                      nCascade, // id
 546                                                                      AliAODVertex::kCascade);
 547     fPrimaryVertex->AddDaughter(vCascade);
 548
 549 //    if (fDebug > 2) {
 550 //      printf("---- Cascade / Cascade Vertex (AOD) : \n");
 551 //      vCascade->Print();
 552 //    }
 553
 554 //    if(esd.GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(const_cast<AliESDEvent*>(&esd), (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production starting form LHC10e without Tender.
 555
 556
 557     // 3 - Add the bachelor track from the cascade
 558
 559     if (!fUsedTrack[idxBachFromCascade]) {
 560
 561       esdCascadeBach->GetPxPyPz(momBach);
 562       esdCascadeBach->GetXYZ(pos);
 563       esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);
 564       esdCascadeBach->GetESDpid(pid);
 565
 566             fUsedTrack[idxBachFromCascade] = kTRUE;
 567             UInt_t selectInfo = 0;
 568             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);
 569             if (fMChandler) fMChandler->SelectParticle(esdCascadeBach->GetLabel());
 570             aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(esdCascadeBach->GetID(),
 571                                                                            esdCascadeBach->GetLabel(),
 572                                                                            momBach,
 573                                                                            kTRUE,
 574                                                                            pos,
 575                                                                            kFALSE, // Why kFALSE for "isDCA" ? FIXME
 576                                                                            covTr,
 577                                                                            (Short_t)esdCascadeBach->GetSign(),
 578                                                                            esdCascadeBach->GetITSClusterMap(),
 579                                                                            pid,
 580                                                                            vCascade,
 581                                                                            kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 582                                                                            vtx->UsesTrack(esdCascadeBach->GetID()),
 583                                                                            AliAODTrack::kSecondary,
 584                                                                            selectInfo);
 585             aodTrack->SetTPCFitMap(esdCascadeBach->GetTPCFitMap());
 586             aodTrack->SetTPCClusterMap(esdCascadeBach->GetTPCClusterMap());
 587             aodTrack->SetTPCSharedMap (esdCascadeBach->GetTPCSharedMap());
 588             aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeBach));
 589             aodTrack->SetTPCPointsF(esdCascadeBach->GetTPCNclsF());
 590             aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeBach->GetTPCCrossedRows()));
 591             fAODTrackRefs->AddAt(aodTrack,idxBachFromCascade);
 592
 593             if (esdCascadeBach->GetSign() > 0) ++fNumberOfPositiveTracks;
 594             aodTrack->ConvertAliPIDtoAODPID();
 595             aodTrack->SetFlags(esdCascadeBach->GetStatus());
 596       SetAODPID(esdCascadeBach,aodTrack,detpid);
 597     }
 598     else {
 599             aodTrack = static_cast<AliAODTrack*>( fAODTrackRefs->At(idxBachFromCascade) );
 600     }
 601
 602     vCascade->AddDaughter(aodTrack);
 603
 604 //    if (fDebug > 4) {
 605 //      printf("---- Cascade / bach dghter : \n");
 606 //      aodTrack->Print();
 607 //    }
 608
 609
 610     // 4 - Add the V0 from the cascade.
 611     // = V0vtx + both pos and neg daughter tracks + the aodV0 itself
 612     //
 613
 614     if ( !fUsedV0[idxV0FromCascade] ) {
 615       // 4.A - if VO structure hasn't been created yet
 616
 617       // 4.A.1 - Create the V0 vertex of the cascade
 618
 619       esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);
 620       esdV0FromCascade->GetPosCov(covVtx);
 621       chi2 = esdV0FromCascade->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3 ?
 622
 623       vV0FromCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex(pos,
 624                                                                          covVtx,
 625                                                                          chi2,
 626                                                                          vCascade,
 627                                                                          idxV0FromCascade, //id of ESDv0
 628                                                                          AliAODVertex::kV0);
 629       // Note:
 630       //    one V0 can be used by several cascades.
 631       // So, one AOD V0 vtx can have several parent vtx.
 632       // This is not directly allowed by AliAODvertex.
 633       // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash
 634       // but to a problem of consistency within AODEvent.
 635       // -> See below paragraph 4.B, for the proposed treatment of such a case.
 636
 637       // Add the vV0FromCascade to the aodVOVtxRefs
 638       fAODV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);
 639
 640
 641       // 4.A.2 - Add the positive tracks from the V0
 642
 643       esdCascadePos->GetPxPyPz(momPos);
 644       esdCascadePos->GetXYZ(pos);
 645       esdCascadePos->GetCovarianceXYZPxPyPz(covTr);
 646       esdCascadePos->GetESDpid(pid);
 647
 648
 649       if (!fUsedTrack[idxPosFromV0Dghter]) {
 650         fUsedTrack[idxPosFromV0Dghter] = kTRUE;
 651
 652         UInt_t selectInfo = 0;
 653         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);
 654         if(fMChandler) fMChandler->SelectParticle(esdCascadePos->GetLabel());
 655         aodTrack = new(tracksArray[fNumberOfTracks++])
 656         AliAODTrack(  esdCascadePos->GetID(),
 657                     esdCascadePos->GetLabel(),
 658                     momPos,
 659                     kTRUE,
 660                     pos,
 661                     kFALSE, // Why kFALSE for "isDCA" ? FIXME
 662                     covTr,
 663                     (Short_t)esdCascadePos->GetSign(),
 664                     esdCascadePos->GetITSClusterMap(),
 665                     pid,
 666                     vV0FromCascade,
 667                     kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 668                     vtx->UsesTrack(esdCascadePos->GetID()),
 669                     AliAODTrack::kSecondary,
 670                     selectInfo);
 671         aodTrack->SetTPCFitMap(esdCascadePos->GetTPCFitMap());
 672         aodTrack->SetTPCClusterMap(esdCascadePos->GetTPCClusterMap());
 673         aodTrack->SetTPCSharedMap (esdCascadePos->GetTPCSharedMap());
 674         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadePos));
 675         aodTrack->SetTPCPointsF(esdCascadePos->GetTPCNclsF());
 676         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadePos->GetTPCCrossedRows()));
 677         fAODTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);
 678
 679         if (esdCascadePos->GetSign() > 0) ++fNumberOfPositiveTracks;
 680         aodTrack->ConvertAliPIDtoAODPID();
 681         aodTrack->SetFlags(esdCascadePos->GetStatus());
 682         SetAODPID(esdCascadePos,aodTrack,detpid);
 683       }
 684       else {
 685         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxPosFromV0Dghter));
 686       }
 687       vV0FromCascade->AddDaughter(aodTrack);
 688
 689
 690       // 4.A.3 - Add the negative tracks from the V0
 691
 692       esdCascadeNeg->GetPxPyPz(momNeg);
 693       esdCascadeNeg->GetXYZ(pos);
 694       esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);
 695       esdCascadeNeg->GetESDpid(pid);
 696
 697
 698       if (!fUsedTrack[idxNegFromV0Dghter]) {
 699         fUsedTrack[idxNegFromV0Dghter] = kTRUE;
 700
 701         UInt_t selectInfo = 0;
 702         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);
 703         if(fMChandler)fMChandler->SelectParticle(esdCascadeNeg->GetLabel());
 704         aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(  esdCascadeNeg->GetID(),
 705                                                       esdCascadeNeg->GetLabel(),
 706                                                       momNeg,
 707                                                       kTRUE,
 708                                                       pos,
 709                                                       kFALSE, // Why kFALSE for "isDCA" ? FIXME
 710                                                       covTr,
 711                                                       (Short_t)esdCascadeNeg->GetSign(),
 712                                                       esdCascadeNeg->GetITSClusterMap(),
 713                                                       pid,
 714                                                       vV0FromCascade,
 715                                                       kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 716                                                       vtx->UsesTrack(esdCascadeNeg->GetID()),
 717                                                       AliAODTrack::kSecondary,
 718                                                       selectInfo);
 719         aodTrack->SetTPCFitMap(esdCascadeNeg->GetTPCFitMap());
 720         aodTrack->SetTPCClusterMap(esdCascadeNeg->GetTPCClusterMap());
 721         aodTrack->SetTPCSharedMap (esdCascadeNeg->GetTPCSharedMap());
 722         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeNeg));
 723         aodTrack->SetTPCPointsF(esdCascadeNeg->GetTPCNclsF());
 724         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeNeg->GetTPCCrossedRows()));
 725         fAODTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);
 726
 727         if (esdCascadeNeg->GetSign() > 0) ++fNumberOfPositiveTracks;
 728         aodTrack->ConvertAliPIDtoAODPID();
 729         aodTrack->SetFlags(esdCascadeNeg->GetStatus());
 730         SetAODPID(esdCascadeNeg,aodTrack,detpid);
 731       }
 732       else {
 733         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxNegFromV0Dghter));
 734       }
 735
 736       vV0FromCascade->AddDaughter(aodTrack);
 737
 738
 739       // 4.A.4 - Add the V0 from cascade to the V0 array
 740
 741       Double_t  dcaV0Daughters      = esdV0FromCascade->GetDcaV0Daughters();
 742       Double_t  dcaV0ToPrimVertex   = esdV0FromCascade->GetD( esd.GetPrimaryVertex()->GetX(),
 743                                                              esd.GetPrimaryVertex()->GetY(),
 744                                                              esd.GetPrimaryVertex()->GetZ() );
 745       esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] );
 746       esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] );
 747
 748       Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
 749       dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD(      esd.GetPrimaryVertex()->GetX(),
 750                                                                   esd.GetPrimaryVertex()->GetY(),
 751                                                                   esd.GetMagneticField())        );
 752       dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD(      esd.GetPrimaryVertex()->GetX(),
 753                                                                   esd.GetPrimaryVertex()->GetY(),
 754                                                                   esd.GetMagneticField())        );
 755
 756       aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0( vV0FromCascade,
 757                                                   dcaV0Daughters,
 758                                                   dcaV0ToPrimVertex,
 759                                                   momPosAtV0vtx,
 760                                                   momNegAtV0vtx,
 761                                                   dcaDaughterToPrimVertex);
 762       // set the aod v0 on-the-fly status
 763       aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());
 764
 765       // Add the aodV0 to the aodVORefs
 766       fAODV0Refs->AddAt(aodV0,idxV0FromCascade);
 767
 768       fUsedV0[idxV0FromCascade] = kTRUE;
 769
 770     } else {
 771       // 4.B - if V0 structure already used
 772
 773       // Note :
 774       //    one V0 can be used by several cascades (frequent in PbPb evts) :
 775       // same V0 which used but attached to different bachelor tracks
 776       // -> aodVORefs and fAODV0VtxRefs are needed.
 777       // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.
 778
 779       vV0FromCascade = static_cast<AliAODVertex*>( fAODV0VtxRefs->At(idxV0FromCascade) );
 780       aodV0          = static_cast<AliAODv0*>    ( fAODV0Refs   ->At(idxV0FromCascade) );
 781
 782       // - Treatment of the parent for such a "re-used" V0 :
 783       // Insert the cascade that reuses the V0 vertex in the lineage chain
 784       // Before : vV0 -> vCascade1 -> vPrimary
 785       //  - Hyp : cascade2 uses the same V0 as cascade1
 786       //  After :  vV0 -> vCascade2 -> vCascade1 -> vPrimary
 787
 788       AliAODVertex *vCascadePreviousParent = static_cast<AliAODVertex*> (vV0FromCascade->GetParent());
 789                         vV0FromCascade->SetParent(vCascade);
 790                         vCascade      ->SetParent(vCascadePreviousParent);
 791
 792 //      if(fDebug > 2)
 793 //        printf("---- Cascade / Lineage insertion\n"
 794 //               "Parent of V0 vtx                 = Cascade vtx %p\n"
 795 //               "Parent of the cascade vtx        = Cascade vtx %p\n"
 796 //               "Parent of the parent cascade vtx = Cascade vtx %p\n",
 797 //               static_cast<void*> (vV0FromCascade->GetParent()),
 798 //               static_cast<void*> (vCascade->GetParent()),
 799 //               static_cast<void*> (vCascadePreviousParent->GetParent()) );
 800
 801     }// end if V0 structure already used
 802
 803 //    if (fDebug > 2) {
 804 //      printf("---- Cascade / V0 vertex: \n");
 805 //      vV0FromCascade->Print();
 806 //    }
 807 //
 808 //    if (fDebug > 4) {
 809 //      printf("---- Cascade / pos dghter : \n");
 810 //                      aodTrack->Print();
 811 //      printf("---- Cascade / neg dghter : \n");
 812 //                      aodTrack->Print();
 813 //      printf("---- Cascade / aodV0 : \n");
 814 //                      aodV0->Print();
 815 //    }
 816
 817     // In any case (used V0 or not), add the V0 vertex to the cascade one.
 818     vCascade->AddDaughter(vV0FromCascade);
 819
 820
 821     // 5 - Add the primary track of the cascade (if any)
 822
 823
 824     // 6 - Add the cascade to the AOD array of cascades
 825
 826     Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd.GetPrimaryVertex()->GetX(),
 827                                                                      esd.GetPrimaryVertex()->GetY(),
 828                                                                      esd.GetMagneticField())        );
 829
 830     Double_t momBachAtCascadeVtx[3]={0.};
 831
 832     esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);
 833
 834     aodCascade = new(cascadesArray[fNumberOfCascades++]) AliAODcascade(  vCascade,
 835                                                           esdCascade->Charge(),
 836                                                           esdCascade->GetDcaXiDaughters(),
 837                                                           -999.,
 838                                                           // DCAXiToPrimVtx -> needs to be calculated   ----|
 839                                                           // doesn't exist at ESD level;
 840                                                           // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)
 841                                                           dcaBachToPrimVertexXY,
 842                                                           momBachAtCascadeVtx,
 843                                                           *aodV0);
 844
 845     if (fDebug > 10) {
 846       printf("---- Cascade / AOD cascade : \n\n");
 847       aodCascade->PrintXi(fPrimaryVertex->GetX(), fPrimaryVertex->GetY(), fPrimaryVertex->GetZ());
 848     }
 849
 850   } // end of the loop on cascades
 851
 852   Cascades().Expand(fNumberOfCascades);
 853 }
 854
 855 //______________________________________________________________________________
 856 void AliAnalysisTaskESDfilter::ConvertV0s(const AliESDEvent& esd)
 857 {
 858   // Access to the AOD container of V0s
 859
 860   AliCodeTimerAuto("",0);
 861
 862   //
 863   // V0s
 864   //
 865
 866   Double_t pos[3] = { 0. };
 867   Double_t chi2(0.0);
 868   Double_t covVtx[6] = { 0. };
 869   Double_t momPos[3]={0.};
 870   Double_t covTr[21]={0.};
 871   Double_t pid[10]={0.};
 872   AliAODTrack* aodTrack(0x0);
 873   AliAODPid* detpid(0x0);
 874   Double_t momNeg[3]={0.};
 875   Double_t momPosAtV0vtx[3]={0.};
 876   Double_t momNegAtV0vtx[3]={0.};
 877
 878   for (Int_t nV0 = 0; nV0 < esd.GetNumberOfV0s(); ++nV0)
 879   {
 880     if (fUsedV0[nV0]) continue; // skip if already added to the AOD
 881
 882     AliESDv0 *v0 = esd.GetV0(nV0);
 883     Int_t posFromV0 = v0->GetPindex();
 884     Int_t negFromV0 = v0->GetNindex();
 885
 886     // V0 selection
 887     //
 888     AliESDVertex *esdVtx   = new AliESDVertex(*(esd.GetPrimaryVertex()));
 889     AliESDtrack  *esdV0Pos = esd.GetTrack(posFromV0);
 890     AliESDtrack  *esdV0Neg = esd.GetTrack(negFromV0);
 891     TList v0objects;
 892     v0objects.AddAt(v0,                      0);
 893     v0objects.AddAt(esdV0Pos,                1);
 894     v0objects.AddAt(esdV0Neg,                2);
 895     v0objects.AddAt(esdVtx,                  3);
 896     UInt_t selectV0 = 0;
 897     if (fV0Filter) {
 898       selectV0 = fV0Filter->IsSelected(&v0objects);
 899       // this is a little awkward but otherwise the
 900       // list wants to access the pointer (delete it)
 901       // again when going out of scope
 902       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 903       esdVtx = 0;
 904       if (!selectV0)
 905         continue;
 906     }
 907     else{
 908       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 909       esdVtx = 0;
 910     }
 911
 912     v0->GetXYZ(pos[0], pos[1], pos[2]);
 913
 914     if (!fOldESDformat) {
 915             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
 916             v0->GetPosCov(covVtx);
 917     } else {
 918             chi2 = -999.;
 919             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 920     }
 921
 922
 923     AliAODVertex * vV0 =
 924           new(Vertices()[fNumberOfVertices++]) AliAODVertex(pos,
 925                                                       covVtx,
 926                                                       chi2,
 927                                                       fPrimaryVertex,
 928                                                       nV0,
 929                                                       AliAODVertex::kV0);
 930     fPrimaryVertex->AddDaughter(vV0);
 931
 932
 933     // Add the positive tracks from the V0
 934
 935
 936     esdV0Pos->GetPxPyPz(momPos);
 937     esdV0Pos->GetXYZ(pos);
 938     esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
 939     esdV0Pos->GetESDpid(pid);
 940
 941     const AliESDVertex *vtx = esd.GetPrimaryVertex();
 942
 943     if (!fUsedTrack[posFromV0]) {
 944             fUsedTrack[posFromV0] = kTRUE;
 945             UInt_t selectInfo = 0;
 946             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);
 947             if(fMChandler)fMChandler->SelectParticle(esdV0Pos->GetLabel());
 948             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Pos->GetID(),
 949                                                     esdV0Pos->GetLabel(),
 950                                                     momPos,
 951                                                     kTRUE,
 952                                                     pos,
 953                                                     kFALSE,
 954                                                     covTr,
 955                                                     (Short_t)esdV0Pos->GetSign(),
 956                                                     esdV0Pos->GetITSClusterMap(),
 957                                                     pid,
 958                                                     vV0,
 959                                                     kTRUE,  // check if this is right
 960                                                     vtx->UsesTrack(esdV0Pos->GetID()),
 961                                                     AliAODTrack::kSecondary,
 962                                                     selectInfo);
 963             aodTrack->SetTPCFitMap(esdV0Pos->GetTPCFitMap());
 964             aodTrack->SetTPCClusterMap(esdV0Pos->GetTPCClusterMap());
 965             aodTrack->SetTPCSharedMap (esdV0Pos->GetTPCSharedMap());
 966             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Pos));
 967             aodTrack->SetTPCPointsF(esdV0Pos->GetTPCNclsF());
 968             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Pos->GetTPCCrossedRows()));
 969             fAODTrackRefs->AddAt(aodTrack,posFromV0);
 970             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
 971             if (esdV0Pos->GetSign() > 0) ++fNumberOfPositiveTracks;
 972             aodTrack->ConvertAliPIDtoAODPID();
 973             aodTrack->SetFlags(esdV0Pos->GetStatus());
 974       SetAODPID(esdV0Pos,aodTrack,detpid);
 975     }
 976     else {
 977             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));
 978             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
 979     }
 980     vV0->AddDaughter(aodTrack);
 981
 982     // Add the negative tracks from the V0
 983
 984     esdV0Neg->GetPxPyPz(momNeg);
 985     esdV0Neg->GetXYZ(pos);
 986     esdV0Neg->GetCovarianceXYZPxPyPz(covTr);
 987     esdV0Neg->GetESDpid(pid);
 988
 989     if (!fUsedTrack[negFromV0]) {
 990             fUsedTrack[negFromV0] = kTRUE;
 991             UInt_t selectInfo = 0;
 992             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);
 993             if(fMChandler)fMChandler->SelectParticle(esdV0Neg->GetLabel());
 994             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Neg->GetID(),
 995                                                     esdV0Neg->GetLabel(),
 996                                                     momNeg,
 997                                                     kTRUE,
 998                                                     pos,
 999                                                     kFALSE,
1000                                                     covTr,
1001                                                     (Short_t)esdV0Neg->GetSign(),
1002                                                     esdV0Neg->GetITSClusterMap(),
1003                                                     pid,
1004                                                     vV0,
1005                                                     kTRUE,  // check if this is right
1006                                                     vtx->UsesTrack(esdV0Neg->GetID()),
1007                                                     AliAODTrack::kSecondary,
1008                                                     selectInfo);
1009             aodTrack->SetTPCFitMap(esdV0Neg->GetTPCFitMap());
1010             aodTrack->SetTPCClusterMap(esdV0Neg->GetTPCClusterMap());
1011             aodTrack->SetTPCSharedMap (esdV0Neg->GetTPCSharedMap());
1012             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Neg));
1013             aodTrack->SetTPCPointsF(esdV0Neg->GetTPCNclsF());
1014             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Neg->GetTPCCrossedRows()));
1015
1016             fAODTrackRefs->AddAt(aodTrack,negFromV0);
1017             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
1018             if (esdV0Neg->GetSign() > 0) ++fNumberOfPositiveTracks;
1019             aodTrack->ConvertAliPIDtoAODPID();
1020             aodTrack->SetFlags(esdV0Neg->GetStatus());
1021       SetAODPID(esdV0Neg,aodTrack,detpid);
1022     }
1023     else {
1024             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(negFromV0));
1025             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
1026     }
1027     vV0->AddDaughter(aodTrack);
1028
1029
1030     // Add the V0 the V0 array as well
1031
1032     Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();
1033     Double_t  dcaV0ToPrimVertex   = v0->GetD(esd.GetPrimaryVertex()->GetX(),
1034                                              esd.GetPrimaryVertex()->GetY(),
1035                                              esd.GetPrimaryVertex()->GetZ());
1036
1037     v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]);
1038     v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]);
1039
1040     Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
1041     dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD(  esd.GetPrimaryVertex()->GetX(),
1042                                                            esd.GetPrimaryVertex()->GetY(),
1043                                                            esd.GetMagneticField()) );
1044     dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD(  esd.GetPrimaryVertex()->GetX(),
1045                                                            esd.GetPrimaryVertex()->GetY(),
1046                                                            esd.GetMagneticField()) );
1047
1048     AliAODv0* aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0(vV0,
1049                                                 dcaV0Daughters,
1050                                                 dcaV0ToPrimVertex,
1051                                                 momPosAtV0vtx,
1052                                                 momNegAtV0vtx,
1053                                                 dcaDaughterToPrimVertex);
1054
1055     // set the aod v0 on-the-fly status
1056     aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
1057   }//End of loop on V0s
1058
1059   V0s().Expand(fNumberOfV0s);
1060 }
1061
1062 //______________________________________________________________________________
1063 void AliAnalysisTaskESDfilter::ConvertTPCOnlyTracks(const AliESDEvent& esd)
1064 {
1065   // Convert TPC only tracks
1066   // Here we have wo hybrid appraoch to remove fakes
1067   // ******* ITSTPC ********
1068   // Uses a cut on the ITS properties to select global tracks
1069   // which are than marked as HybdridITSTPC for the remainder
1070   // the TPC only tracks are flagged as HybridITSTPConly.
1071   // Note, in order not to get fakes back in the TPC cuts, one needs
1072   // two "ITS" cuts one tight (1) (to throw out fakes) and one lose (2) (to NOT flag the trakcs in the TPC only)
1073   // using cut number (3)
1074   // so fHybridFilterMask == (1)|(2) fTPCFilterMask = (3), Usercode needs to slect with mask = (1)|(3) and track->IsHybridITSTPC()
1075   // ******* TPC ********
1076   // Here, only TPC tracks are flagged that pass the tight ITS cuts and tracks that pass the TPC cuts and NOT the loose ITS cuts
1077   // the ITS cuts neeed to be added to the filter as extra cuts, since here the selections info is reset in the global and put to the TPC only track
1078
1079   AliCodeTimerAuto("",0);
1080
1081   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1082   for(int it = 0;it < fNumberOfTracks;++it)
1083   {
1084     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1085     if(!tr)continue;
1086     UInt_t map = tr->GetFilterMap();
1087     if(map&fTPCConstrainedFilterMask){
1088       // we only reset the track select ionfo, no deletion...
1089       tr->SetFilterMap(map&~fTPCConstrainedFilterMask);
1090     }
1091     if(map&fHybridFilterMaskTPCCG){
1092       // this is one part of the hybrid tracks
1093       // the others not passing the selection will be TPC only selected below
1094       tr->SetIsHybridTPCConstrainedGlobal(kTRUE);
1095     }
1096   }
1097   // Loop over the ESD trcks and pick out the tracks passing TPC only cuts
1098
1099
1100   const AliESDVertex *vtxSPD = esd.GetPrimaryVertexSPD();
1101   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1102
1103   Double_t pos[3] = { 0. };
1104   Double_t covTr[21]={0.};
1105   Double_t pid[10]={0.};
1106
1107
1108   Double_t p[3] = { 0. };
1109
1110   Double_t pDCA[3] = { 0. }; // momentum at DCA
1111   Double_t rDCA[3] = { 0. }; // position at DCA
1112   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1113   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1114
1115
1116   AliAODTrack* aodTrack(0x0);
1117   //  AliAODPid* detpid(0x0);
1118
1119   // account for change in pT after the constraint
1120   Float_t ptMax = 1E10;
1121   Float_t ptMin = 0;
1122   for(int i = 0;i<32;i++){
1123     if(fTPCConstrainedFilterMask&(1<<i)){
1124       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1125       Float_t tmp1= 0,tmp2 = 0;
1126       cuts->GetPtRange(tmp1,tmp2);
1127       if(tmp1>ptMin)ptMin=tmp1;
1128       if(tmp2<ptMax)ptMax=tmp2;
1129     }
1130   }
1131
1132   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1133   {
1134     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1135
1136     UInt_t selectInfo = 0;
1137     Bool_t isHybridITSTPC = false;
1138     //
1139     // Track selection
1140     if (fTrackFilter) {
1141       selectInfo = fTrackFilter->IsSelected(esdTrack);
1142     }
1143
1144     if(!(selectInfo&fHybridFilterMaskTPCCG)){
1145       // not already selected tracks, use second part of hybrid tracks
1146       isHybridITSTPC = true;
1147       // too save space one could only store these...
1148     }
1149
1150     selectInfo &= fTPCConstrainedFilterMask;
1151     if (!selectInfo)continue;
1152     if (fWriteHybridTPCCOnly&&!isHybridITSTPC)continue; // write only complementary tracks
1153     // create a tpc only tracl
1154     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(&esd),esdTrack->GetID());
1155     if(!track) continue;
1156
1157     if(track->Pt()>0.)
1158     {
1159       // only constrain tracks above threshold
1160       AliExternalTrackParam exParam;
1161       // take the B-field from the ESD, no 3D fieldMap available at this point
1162       Bool_t relate = false;
1163       relate = track->RelateToVertexTPC(vtxSPD,esd.GetMagneticField(),kVeryBig,&exParam);
1164       if(!relate){
1165         delete track;
1166         continue;
1167       }
1168       // fetch the track parameters at the DCA (unconstraint)
1169       if(track->GetTPCInnerParam()){
1170         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1171         track->GetTPCInnerParam()->GetXYZ(rDCA);
1172       }
1173       // get the DCA to the vertex:
1174       track->GetImpactParametersTPC(dDCA,cDCA);
1175       // set the constrained parameters to the track
1176       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1177     }
1178
1179     track->GetPxPyPz(p);
1180
1181     Float_t pT = track->Pt();
1182     if(pT<ptMin||pT>ptMax){
1183       delete track;
1184       continue;
1185     }
1186
1187     //
1188
1189
1190     track->GetXYZ(pos);
1191     track->GetCovarianceXYZPxPyPz(covTr);
1192     esdTrack->GetESDpid(pid);// original PID
1193
1194     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1195     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((track->GetID()+1)*-1,
1196                                                             track->GetLabel(),
1197                                                             p,
1198                                                             kTRUE,
1199                                                             pos,
1200                                                             kFALSE,
1201                                                             covTr,
1202                                                             (Short_t)track->GetSign(),
1203                                                             track->GetITSClusterMap(),
1204                                                             pid,
1205                                                             fPrimaryVertex,
1206                                                             kTRUE, // check if this is right
1207                                                             vtx->UsesTrack(track->GetID()),
1208                                                             AliAODTrack::kPrimary,
1209                                                             selectInfo);
1210     aodTrack->SetIsHybridTPCConstrainedGlobal(isHybridITSTPC);
1211     aodTrack->SetTPCFitMap(track->GetTPCFitMap());
1212     aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
1213     aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
1214     aodTrack->SetIsTPCConstrained(kTRUE);
1215     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack)); // original track
1216     // set the DCA values to the AOD track
1217     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1218     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1219     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1220
1221     aodTrack->SetFlags(track->GetStatus());
1222     aodTrack->SetTPCPointsF(track->GetTPCNclsF());
1223     aodTrack->SetTPCNCrossedRows(UShort_t(track->GetTPCCrossedRows()));
1224
1225     //Perform progagation of tracks if needed
1226     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1227     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1228
1229     // do not duplicate PID information
1230     //    aodTrack->ConvertAliPIDtoAODPID();
1231     //    SetAODPID(esdTrack,aodTrack,detpid);
1232
1233     delete track;
1234   } // end of loop on tracks
1235
1236 }
1237
1238 //______________________________________________________________________________
1239 void AliAnalysisTaskESDfilter::ConvertGlobalConstrainedTracks(const AliESDEvent& esd)
1240 {
1241
1242   // Here we have the option to store the complement from global constraint information
1243   // to tracks passing tight cuts (1) in order not to get fakes back in, one needs
1244   // two sets of cuts one tight (1) (to throw out fakes) and one lose (2) (fakes/bad tracks would pass (2) but not (1))
1245   // using cut number (3) selects the tracks that complement (1) e.g. tracks witout ITS refit or cluster requirement
1246
1247
1248   AliCodeTimerAuto("",0);
1249
1250   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1251   for(int it = 0;it < fNumberOfTracks;++it)
1252   {
1253     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1254     if(!tr)continue;
1255     UInt_t map = tr->GetFilterMap();
1256     if(map&fGlobalConstrainedFilterMask){
1257       // we only reset the track select info, no deletion...
1258       // mask reset mask in case track is already taken
1259       tr->SetFilterMap(map&~fGlobalConstrainedFilterMask);
1260     }
1261     if(map&fHybridFilterMaskGCG){
1262       // this is one part of the hybrid tracks
1263       // the others not passing the selection will be the ones selected below
1264       tr->SetIsHybridGlobalConstrainedGlobal(kTRUE);
1265     }
1266   }
1267   // Loop over the ESD trcks and pick out the tracks passing the GlobalConstraint cuts
1268
1269
1270   Double_t pos[3] = { 0. };
1271   Double_t covTr[21]={0.};
1272   Double_t pid[10]={0.};
1273   Double_t p[3] = { 0. };
1274
1275   Double_t pDCA[3] = { 0. }; // momentum at DCA
1276   Double_t rDCA[3] = { 0. }; // position at DCA
1277   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1278   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1279
1280
1281   AliAODTrack* aodTrack(0x0);
1282   AliAODPid* detpid(0x0);
1283   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1284
1285   // account for change in pT after the constraint
1286   Float_t ptMax = 1E10;
1287   Float_t ptMin = 0;
1288   for(int i = 0;i<32;i++){
1289     if(fGlobalConstrainedFilterMask&(1<<i)){
1290       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1291       Float_t tmp1= 0,tmp2 = 0;
1292       cuts->GetPtRange(tmp1,tmp2);
1293       if(tmp1>ptMin)ptMin=tmp1;
1294       if(tmp2<ptMax)ptMax=tmp2;
1295     }
1296   }
1297
1298
1299
1300  for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1301   {
1302     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1303     const AliExternalTrackParam * exParamGC = esdTrack->GetConstrainedParam();
1304     if(!exParamGC)continue;
1305
1306     UInt_t selectInfo = 0;
1307     Bool_t isHybridGC = false;
1308
1309     //
1310     // Track selection
1311     if (fTrackFilter) {
1312       selectInfo = fTrackFilter->IsSelected(esdTrack);
1313     }
1314
1315
1316     if(!(selectInfo&fHybridFilterMaskGCG))isHybridGC = true;
1317     if (fWriteHybridGCOnly&&!isHybridGC)continue; // write only complementary tracks
1318
1319     selectInfo &= fGlobalConstrainedFilterMask;
1320     if (!selectInfo)continue;
1321     // fetch the track parameters at the DCA (unconstrained)
1322     esdTrack->GetPxPyPz(pDCA);
1323     esdTrack->GetXYZ(rDCA);
1324     // get the DCA to the vertex:
1325     esdTrack->GetImpactParameters(dDCA,cDCA);
1326
1327     if (!esdTrack->GetConstrainedPxPyPz(p)) continue;
1328
1329
1330     Float_t pT = exParamGC->Pt();
1331     if(pT<ptMin||pT>ptMax){
1332       continue;
1333     }
1334
1335
1336     esdTrack->GetConstrainedXYZ(pos);
1337     exParamGC->GetCovarianceXYZPxPyPz(covTr);
1338     esdTrack->GetESDpid(pid);
1339     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1340     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((esdTrack->GetID()+1)*-1,
1341                                                             esdTrack->GetLabel(),
1342                                                             p,
1343                                                             kTRUE,
1344                                                             pos,
1345                                                             kFALSE,
1346                                                             covTr,
1347                                                             (Short_t)esdTrack->GetSign(),
1348                                                             esdTrack->GetITSClusterMap(),
1349                                                             pid,
1350                                                             fPrimaryVertex,
1351                                                             kTRUE, // check if this is right
1352                                                             vtx->UsesTrack(esdTrack->GetID()),
1353                                                             AliAODTrack::kPrimary,
1354                                                             selectInfo);
1355     aodTrack->SetIsHybridGlobalConstrainedGlobal(isHybridGC);
1356     aodTrack->SetIsGlobalConstrained(kTRUE);
1357     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1358     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1359     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1360     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1361
1362
1363     // set the DCA values to the AOD track
1364     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1365     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1366     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1367
1368     aodTrack->SetFlags(esdTrack->GetStatus());
1369     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1370     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1371
1372     if(isHybridGC){
1373       // only copy AOD information for hybrid, no duplicate information
1374       aodTrack->ConvertAliPIDtoAODPID();
1375       SetAODPID(esdTrack,aodTrack,detpid);
1376     }
1377
1378     //Perform progagation of tracks if needed
1379     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1380     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1381   } // end of loop on tracks
1382
1383 }
1384
1385
1386 //______________________________________________________________________________
1387 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)
1388 {
1389   // Tracks (primary and orphan)
1390
1391   AliCodeTimerAuto("",0);
1392
1393   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));
1394
1395   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1396   Double_t p[3] = { 0. };
1397   Double_t pos[3] = { 0. };
1398   Double_t covTr[21] = { 0. };
1399   Double_t pid[10] = { 0. };
1400   AliAODTrack* aodTrack(0x0);
1401   AliAODPid* detpid(0x0);
1402
1403   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1404   {
1405     if (fUsedTrack[nTrack]) continue;
1406
1407     AliESDtrack *esdTrack = esd.GetTrack(nTrack);
1408     UInt_t selectInfo = 0;
1409     //
1410     // Track selection
1411     if (fTrackFilter) {
1412             selectInfo = fTrackFilter->IsSelected(esdTrack);
1413             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;
1414     }
1415
1416
1417     esdTrack->GetPxPyPz(p);
1418     esdTrack->GetXYZ(pos);
1419     esdTrack->GetCovarianceXYZPxPyPz(covTr);
1420     esdTrack->GetESDpid(pid);
1421     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1422     fPrimaryVertex->AddDaughter(aodTrack =
1423                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),
1424                                                             esdTrack->GetLabel(),
1425                                                             p,
1426                                                             kTRUE,
1427                                                             pos,
1428                                                             kFALSE,
1429                                                             covTr,
1430                                                             (Short_t)esdTrack->GetSign(),
1431                                                             esdTrack->GetITSClusterMap(),
1432                                                             pid,
1433                                                             fPrimaryVertex,
1434                                                             kTRUE, // check if this is right
1435                                                             vtx->UsesTrack(esdTrack->GetID()),
1436                                                             AliAODTrack::kPrimary,
1437                                                             selectInfo)
1438                          );
1439     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1440     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1441     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1442     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1443     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1444     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1445     if(esdTrack->IsEMCAL()) aodTrack->SetEMCALcluster(esdTrack->GetEMCALcluster());
1446     if(esdTrack->IsPHOS())  aodTrack->SetPHOScluster(esdTrack->GetPHOScluster());
1447
1448     //Perform progagation of tracks if needed
1449     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1450     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1451
1452     fAODTrackRefs->AddAt(aodTrack, nTrack);
1453
1454
1455     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;
1456     aodTrack->SetFlags(esdTrack->GetStatus());
1457     aodTrack->ConvertAliPIDtoAODPID();
1458     SetAODPID(esdTrack,aodTrack,detpid);
1459   } // end of loop on tracks
1460 }
1461
1462 //______________________________________________________________________________
1463 void AliAnalysisTaskESDfilter::PropagateTrackToEMCal(AliESDtrack *esdTrack)
1464 {
1465   Double_t trkPos[3] = {0.,0.,0.};
1466   Double_t EMCalEta=-999, EMCalPhi=-999;
1467   Double_t trkphi = esdTrack->Phi()*TMath::RadToDeg();
1468   if(TMath::Abs(esdTrack->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
1469     {
1470       AliExternalTrackParam *trkParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1471       if(trkParam)
1472         {
1473           AliExternalTrackParam trkParamTmp(*trkParam);
1474           if(AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, 430, esdTrack->GetMass(), 20, kTRUE, 0.8, -1))
1475             {
1476               trkParamTmp.GetXYZ(trkPos);
1477               TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1478               EMCalEta = trkPosVec.Eta();
1479               EMCalPhi = trkPosVec.Phi();
1480               if(EMCalPhi<0)  EMCalPhi += 2*TMath::Pi();
1481               esdTrack->SetTrackPhiEtaOnEMCal(EMCalPhi,EMCalEta);
1482             }
1483         }
1484     }
1485 }
1486
1487 //______________________________________________________________________________
1488 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)
1489 {
1490 // Convert PMD Clusters
1491   AliCodeTimerAuto("",0);
1492   Int_t jPmdClusters=0;
1493   // Access to the AOD container of PMD clusters
1494   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
1495   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {
1496     // file pmd clusters, to be revised!
1497     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);
1498     Int_t nLabel = 0;
1499     Int_t *label = 0x0;
1500     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
1501     Double_t pidPmd[13] = { 0.}; // to be revised!
1502     // type not set!
1503     // assoc cluster not set
1504     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
1505   }
1506 }
1507
1508
1509 //______________________________________________________________________________
1510 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)
1511 {
1512 // Convert Calorimeter Clusters
1513   AliCodeTimerAuto("",0);
1514
1515   // Access to the AOD container of clusters
1516   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
1517   Int_t jClusters(0);
1518
1519   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
1520
1521     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
1522
1523     Int_t  id        = cluster->GetID();
1524     Int_t  nLabel    = cluster->GetNLabels();
1525     Int_t *labels    = cluster->GetLabels();
1526     if(labels){
1527                   for(int i = 0;i < nLabel;++i){
1528                           if(fMChandler)fMChandler->SelectParticle(labels[i]);
1529                   }
1530           }
1531
1532     Float_t energy = cluster->E();
1533     Float_t posF[3] = { 0.};
1534     cluster->GetPosition(posF);
1535
1536     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
1537                                                                                       nLabel,
1538                                                                                       labels,
1539                                                                                       energy,
1540                                                                                       posF,
1541                                                                                       NULL,
1542                                                                                       cluster->GetType(),0);
1543
1544     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
1545                                 cluster->GetDispersion(),
1546                                 cluster->GetM20(), cluster->GetM02(),
1547                                 cluster->GetEmcCpvDistance(),
1548                                 cluster->GetNExMax(),cluster->GetTOF()) ;
1549
1550     caloCluster->SetPIDFromESD(cluster->GetPID());
1551     caloCluster->SetNCells(cluster->GetNCells());
1552     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
1553     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
1554
1555     caloCluster->SetTrackDistance(cluster->GetTrackDx(), cluster->GetTrackDz());
1556
1557     Int_t nMatchCount = 0;
1558     TArrayI* matchedT =         cluster->GetTracksMatched();
1559     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
1560       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
1561         Int_t iESDtrack = matchedT->At(im);;
1562         if (fAODTrackRefs->At(iESDtrack) != 0) {
1563           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
1564           nMatchCount++;
1565         }
1566       }
1567     }
1568     if(nMatchCount==0)
1569       caloCluster->SetTrackDistance(-999,-999);
1570
1571   }
1572   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
1573 }
1574
1575 //______________________________________________________________________________
1576 void AliAnalysisTaskESDfilter::ConvertCaloTrigger(TString calo, const AliESDEvent& esd)
1577 {
1578         AliCodeTimerAuto("",0);
1579
1580         if (calo == "PHOS")
1581         {
1582           AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1583           AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1584
1585           aodTrigger.Allocate(esdTrigger.GetEntries());
1586           esdTrigger.Reset();
1587
1588           Float_t a;
1589           Int_t tmod,tabsId;
1590
1591           while (esdTrigger.Next()) {
1592             esdTrigger.GetPosition(tmod,tabsId);
1593             esdTrigger.GetAmplitude(a);
1594             aodTrigger.Add(tmod,tabsId,a,0.,(Int_t*)NULL,0,0,0);
1595           }
1596
1597           return;
1598         }
1599
1600         AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1601
1602         if (aodHandler)
1603         {
1604                 TTree *aodTree = aodHandler->GetTree();
1605
1606                 if (aodTree)
1607                 {
1608                         Int_t *type = esd.GetCaloTriggerType();
1609
1610                         for (Int_t i = 0; i < 15; i++)
1611                         {
1612                                 aodTree->GetUserInfo()->Add(new TParameter<int>(Form("EMCALCaloTrigger%d",i), type[i]));
1613                         }
1614                 }
1615         }
1616
1617         AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1618
1619         AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1620
1621         aodTrigger.Allocate(esdTrigger.GetEntries());
1622
1623         esdTrigger.Reset();
1624         while (esdTrigger.Next())
1625         {
1626                 Int_t px, py, ts, nTimes, times[10], b;
1627                 Float_t a, t;
1628
1629                 esdTrigger.GetPosition(px, py);
1630
1631                 esdTrigger.GetAmplitude(a);
1632                 esdTrigger.GetTime(t);
1633
1634                 esdTrigger.GetL0Times(times);
1635                 esdTrigger.GetNL0Times(nTimes);
1636
1637                 esdTrigger.GetL1TimeSum(ts);
1638
1639                 esdTrigger.GetTriggerBits(b);
1640
1641                 aodTrigger.Add(px, py, a, t, times, nTimes, ts, b);
1642         }
1643
1644         for (int i = 0; i < 4; i++) aodTrigger.SetL1Threshold(i, esdTrigger.GetL1Threshold(i));
1645
1646         Int_t v0[2] =
1647         {
1648                 esdTrigger.GetL1V0(0),
1649                 esdTrigger.GetL1V0(1)
1650         };
1651
1652         aodTrigger.SetL1V0(v0);
1653         aodTrigger.SetL1FrameMask(esdTrigger.GetL1FrameMask());
1654 }
1655
1656 //______________________________________________________________________________
1657 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)
1658 {
1659 // Convert EMCAL Cells
1660   AliCodeTimerAuto("",0);
1661   // fill EMCAL cell info
1662   if (esd.GetEMCALCells()) { // protection against missing ESD information
1663     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
1664     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
1665
1666     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1667     aodEMcells.CreateContainer(nEMcell);
1668     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
1669     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
1670       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),
1671                          esdEMcells.GetTime(iCell), esdEMcells.GetMCLabel(iCell), esdEMcells.GetEFraction(iCell));
1672     }
1673     aodEMcells.Sort();
1674   }
1675 }
1676
1677 //______________________________________________________________________________
1678 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)
1679 {
1680 // Convert PHOS Cells
1681   AliCodeTimerAuto("",0);
1682   // fill PHOS cell info
1683   if (esd.GetPHOSCells()) { // protection against missing ESD information
1684     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
1685     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
1686
1687     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
1688     aodPHcells.CreateContainer(nPHcell);
1689     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
1690     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
1691       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),
1692                          esdPHcells.GetTime(iCell), esdPHcells.GetMCLabel(iCell), esdPHcells.GetEFraction(iCell));
1693     }
1694     aodPHcells.Sort();
1695   }
1696 }
1697
1698 //______________________________________________________________________________
1699 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)
1700 {
1701   // tracklets
1702   AliCodeTimerAuto("",0);
1703
1704   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());
1705   const AliMultiplicity *mult = esd.GetMultiplicity();
1706   if (mult) {
1707     if (mult->GetNumberOfTracklets()>0) {
1708       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
1709
1710       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
1711         if(fMChandler){
1712           fMChandler->SelectParticle(mult->GetLabel(n, 0));
1713           fMChandler->SelectParticle(mult->GetLabel(n, 1));
1714         }
1715         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1716       }
1717     }
1718   } else {
1719     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
1720   }
1721 }
1722
1723 //______________________________________________________________________________
1724 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)
1725 {
1726   AliCodeTimerAuto("",0);
1727
1728   // Kinks: it is a big mess the access to the information in the kinks
1729   // The loop is on the tracks in order to find the mother and daugther of each kink
1730
1731   Double_t covTr[21]={0.};
1732   Double_t pid[10]={0.};
1733   AliAODPid* detpid(0x0);
1734
1735   fNumberOfKinks = esd.GetNumberOfKinks();
1736
1737   const AliESDVertex* vtx = esd.GetPrimaryVertex();
1738
1739   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
1740   {
1741     AliESDtrack * esdTrack = esd.GetTrack(iTrack);
1742
1743     Int_t ikink = esdTrack->GetKinkIndex(0);
1744
1745     if (ikink && fNumberOfKinks) {
1746             // Negative kink index: mother, positive: daughter
1747
1748             // Search for the second track of the kink
1749
1750             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {
1751
1752         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);
1753
1754         Int_t jkink = esdTrack1->GetKinkIndex(0);
1755
1756         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
1757
1758           // The two tracks are from the same kink
1759
1760           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
1761
1762           Int_t imother = -1;
1763           Int_t idaughter = -1;
1764
1765           if (ikink<0 && jkink>0) {
1766
1767             imother = iTrack;
1768             idaughter = jTrack;
1769           }
1770           else if (ikink>0 && jkink<0) {
1771
1772             imother = jTrack;
1773             idaughter = iTrack;
1774           }
1775           else {
1776             //                  cerr << "Error: Wrong combination of kink indexes: "
1777             //                       << ikink << " " << jkink << endl;
1778             continue;
1779           }
1780
1781           // Add the mother track if it passed primary track selection cuts
1782
1783           AliAODTrack * mother = NULL;
1784
1785           UInt_t selectInfo = 0;
1786           if (fTrackFilter) {
1787             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));
1788             if (!selectInfo) continue;
1789           }
1790
1791           if (!fUsedTrack[imother]) {
1792
1793             fUsedTrack[imother] = kTRUE;
1794
1795             AliESDtrack *esdTrackM = esd.GetTrack(imother);
1796             Double_t p[3] = { 0. };
1797             Double_t pos[3] = { 0. };
1798             esdTrackM->GetPxPyPz(p);
1799             esdTrackM->GetXYZ(pos);
1800             esdTrackM->GetCovarianceXYZPxPyPz(covTr);
1801             esdTrackM->GetESDpid(pid);
1802             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());
1803             mother =
1804             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),
1805                                                esdTrackM->GetLabel(),
1806                                                p,
1807                                                kTRUE,
1808                                                pos,
1809                                                kFALSE,
1810                                                covTr,
1811                                                (Short_t)esdTrackM->GetSign(),
1812                                                esdTrackM->GetITSClusterMap(),
1813                                                pid,
1814                                                fPrimaryVertex,
1815                                                kTRUE, // check if this is right
1816                                                vtx->UsesTrack(esdTrack->GetID()),
1817                                                AliAODTrack::kPrimary,
1818                                                selectInfo);
1819             mother->SetTPCFitMap(esdTrackM->GetTPCFitMap());
1820             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());
1821             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());
1822             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));
1823             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());
1824             mother->SetTPCNCrossedRows(UShort_t(esdTrackM->GetTPCCrossedRows()));
1825
1826             fAODTrackRefs->AddAt(mother, imother);
1827
1828             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;
1829             mother->SetFlags(esdTrackM->GetStatus());
1830             mother->ConvertAliPIDtoAODPID();
1831             fPrimaryVertex->AddDaughter(mother);
1832             mother->ConvertAliPIDtoAODPID();
1833             SetAODPID(esdTrackM,mother,detpid);
1834           }
1835           else {
1836             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1837             //                       << " track " << imother << " has already been used!" << endl;
1838           }
1839
1840           // Add the kink vertex
1841           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);
1842
1843           AliAODVertex * vkink =
1844           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),
1845                                                   NULL,
1846                                                   0.,
1847                                                   mother,
1848                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!
1849                                                   AliAODVertex::kKink);
1850           // Add the daughter track
1851
1852           AliAODTrack * daughter = NULL;
1853
1854           if (!fUsedTrack[idaughter]) {
1855
1856             fUsedTrack[idaughter] = kTRUE;
1857
1858             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);
1859             Double_t p[3] = { 0. };
1860             Double_t pos[3] = { 0. };
1861
1862             esdTrackD->GetPxPyPz(p);
1863             esdTrackD->GetXYZ(pos);
1864             esdTrackD->GetCovarianceXYZPxPyPz(covTr);
1865             esdTrackD->GetESDpid(pid);
1866             selectInfo = 0;
1867             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);
1868             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());
1869             daughter =
1870             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),
1871                                                esdTrackD->GetLabel(),
1872                                                p,
1873                                                kTRUE,
1874                                                pos,
1875                                                kFALSE,
1876                                                covTr,
1877                                                (Short_t)esdTrackD->GetSign(),
1878                                                esdTrackD->GetITSClusterMap(),
1879                                                pid,
1880                                                vkink,
1881                                                kTRUE, // check if this is right
1882                                                vtx->UsesTrack(esdTrack->GetID()),
1883                                                AliAODTrack::kSecondary,
1884                                                selectInfo);
1885             daughter->SetTPCFitMap(esdTrackD->GetTPCFitMap());
1886             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());
1887             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());
1888             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());
1889             daughter->SetTPCNCrossedRows(UShort_t(esdTrackD->GetTPCCrossedRows()));
1890             fAODTrackRefs->AddAt(daughter, idaughter);
1891
1892             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;
1893             daughter->SetFlags(esdTrackD->GetStatus());
1894             daughter->ConvertAliPIDtoAODPID();
1895             vkink->AddDaughter(daughter);
1896             daughter->ConvertAliPIDtoAODPID();
1897             SetAODPID(esdTrackD,daughter,detpid);
1898           }
1899           else {
1900             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1901             //                       << " track " << idaughter << " has already been used!" << endl;
1902           }
1903         }
1904             }
1905     }
1906   }
1907 }
1908
1909 //______________________________________________________________________________
1910 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)
1911 {
1912   AliCodeTimerAuto("",0);
1913
1914   // Access to the AOD container of vertices
1915   fNumberOfVertices = 0;
1916
1917   Double_t pos[3] = { 0. };
1918   Double_t covVtx[6] = { 0. };
1919
1920   // Add primary vertex. The primary tracks will be defined
1921   // after the loops on the composite objects (V0, cascades, kinks)
1922   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1923
1924   vtx->GetXYZ(pos); // position
1925   vtx->GetCovMatrix(covVtx); //covariance matrix
1926
1927   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])
1928   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1929   fPrimaryVertex->SetName(vtx->GetName());
1930   fPrimaryVertex->SetTitle(vtx->GetTitle());
1931   fPrimaryVertex->SetBC(vtx->GetBC());
1932
1933   TString vtitle = vtx->GetTitle();
1934   if (!vtitle.Contains("VertexerTracks"))
1935     fPrimaryVertex->SetNContributors(vtx->GetNContributors());
1936
1937   if (fDebug > 0) fPrimaryVertex->Print();
1938
1939   // Add SPD "main" vertex
1940   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
1941   vtxS->GetXYZ(pos); // position
1942   vtxS->GetCovMatrix(covVtx); //covariance matrix
1943   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])
1944   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
1945   mVSPD->SetName(vtxS->GetName());
1946   mVSPD->SetTitle(vtxS->GetTitle());
1947   mVSPD->SetNContributors(vtxS->GetNContributors());
1948
1949   // Add SPD pileup vertices
1950   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
1951   {
1952     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
1953     vtxP->GetXYZ(pos); // position
1954     vtxP->GetCovMatrix(covVtx); //covariance matrix
1955     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])
1956     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
1957     pVSPD->SetName(vtxP->GetName());
1958     pVSPD->SetTitle(vtxP->GetTitle());
1959     pVSPD->SetNContributors(vtxP->GetNContributors());
1960     pVSPD->SetBC(vtxP->GetBC());
1961   }
1962
1963   // Add TRK pileup vertices
1964   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
1965   {
1966     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
1967     vtxP->GetXYZ(pos); // position
1968     vtxP->GetCovMatrix(covVtx); //covariance matrix
1969     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])
1970     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
1971     pVTRK->SetName(vtxP->GetName());
1972     pVTRK->SetTitle(vtxP->GetTitle());
1973     pVTRK->SetNContributors(vtxP->GetNContributors());
1974     pVTRK->SetBC(vtxP->GetBC());
1975   }
1976 }
1977
1978 //______________________________________________________________________________
1979 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)
1980 {
1981   // Convert VZERO data
1982   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
1983   *vzeroData = *(esd.GetVZEROData());
1984 }
1985
1986 //______________________________________________________________________________
1987 void AliAnalysisTaskESDfilter::ConvertTZERO(const AliESDEvent& esd)
1988 {
1989   // Convert TZERO data
1990   const AliESDTZERO* esdTzero = esd.GetESDTZERO();
1991   AliAODTZERO* aodTzero = AODEvent()->GetTZEROData();
1992
1993   for (Int_t icase=0; icase<3; icase++){
1994     aodTzero->SetT0TOF(    icase, esdTzero->GetT0TOF(icase));
1995     aodTzero->SetT0TOFbest(icase, esdTzero->GetT0TOFbest(icase));
1996   }
1997   aodTzero->SetBackgroundFlag(esdTzero->GetBackgroundFlag());
1998   aodTzero->SetPileupFlag(esdTzero->GetPileupFlag());
1999   aodTzero->SetSatelliteFlag(esdTzero->GetSatellite());
2000
2001   Float_t rawTime[24];
2002   for(Int_t ipmt=0; ipmt<24; ipmt++)
2003     rawTime[ipmt] = esdTzero->GetTimeFull(ipmt,0);
2004
2005   Int_t idxOfFirstPmtA = -1,       idxOfFirstPmtC = -1;
2006   Float_t timeOfFirstPmtA = 9999, timeOfFirstPmtC = 9999;
2007   for(int ipmt=0;  ipmt<12; ipmt++){
2008     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtC && rawTime[ipmt]!=0){
2009       timeOfFirstPmtC = rawTime[ipmt];
2010       idxOfFirstPmtC  = ipmt;
2011     }
2012   }
2013   for(int ipmt=12; ipmt<24; ipmt++){
2014     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtA && rawTime[ipmt]!=0 ){
2015       timeOfFirstPmtA = rawTime[ipmt];
2016       idxOfFirstPmtA  = ipmt;
2017     }
2018   }
2019
2020   if(idxOfFirstPmtA != -1 && idxOfFirstPmtC != -1){
2021     //speed of light in cm/ns   TMath::C()*1e-7
2022     Float_t vertexraw = TMath::C()*1e-7 * (rawTime[idxOfFirstPmtA] - rawTime[idxOfFirstPmtC])/2;
2023     aodTzero->SetT0VertexRaw( vertexraw );
2024   }else{
2025     aodTzero->SetT0VertexRaw(99999);
2026   }
2027
2028   aodTzero->SetT0zVertex(esdTzero->GetT0zVertex());
2029 }
2030
2031
2032 //______________________________________________________________________________
2033 void AliAnalysisTaskESDfilter::ConvertZDC(const AliESDEvent& esd)
2034 {
2035   // Convert ZDC data
2036   AliESDZDC* esdZDC = esd.GetZDCData();
2037
2038   const Double_t zem1Energy = esdZDC->GetZEM1Energy();
2039   const Double_t zem2Energy = esdZDC->GetZEM2Energy();
2040
2041   const Double_t *towZNC = esdZDC->GetZNCTowerEnergy();
2042   const Double_t *towZPC = esdZDC->GetZPCTowerEnergy();
2043   const Double_t *towZNA = esdZDC->GetZNATowerEnergy();
2044   const Double_t *towZPA = esdZDC->GetZPATowerEnergy();
2045   const Double_t *towZNCLG = esdZDC->GetZNCTowerEnergyLR();
2046   const Double_t *towZNALG = esdZDC->GetZNATowerEnergyLR();
2047
2048   AliAODZDC* zdcAOD = AODEvent()->GetZDCData();
2049
2050   zdcAOD->SetZEM1Energy(zem1Energy);
2051   zdcAOD->SetZEM2Energy(zem2Energy);
2052   zdcAOD->SetZNCTowers(towZNC, towZNCLG);
2053   zdcAOD->SetZNATowers(towZNA, towZNALG);
2054   zdcAOD->SetZPCTowers(towZPC);
2055   zdcAOD->SetZPATowers(towZPA);
2056
2057   zdcAOD->SetZDCParticipants(esdZDC->GetZDCParticipants(), esdZDC->GetZDCPartSideA(), esdZDC->GetZDCPartSideC());
2058   zdcAOD->SetZDCImpactParameter(esdZDC->GetImpactParameter(), esdZDC->GetImpactParamSideA(),
2059         esdZDC->GetImpactParamSideC());
2060   zdcAOD->SetZDCTDCSum(esdZDC->GetZNTDCSum(0));
2061   zdcAOD->SetZDCTDCDiff(esdZDC->GetZNTDCDiff(0));
2062
2063 }
2064
2065 //_______________________________________________________________________________________________________________________________________
2066 Int_t AliAnalysisTaskESDfilter::ConvertHMPID(const AliESDEvent& esd) // clm
2067 {
2068   //
2069   // Convtert ESD HMPID info to AOD and return the number of good tracks with HMPID signal.
2070   // We need to return an int since there is no signal counter in the ESD.
2071   //
2072
2073   AliCodeTimerAuto("",0);
2074
2075   Int_t cntHmpidGoodTracks = 0;
2076
2077   Float_t  xMip = 0;
2078   Float_t  yMip = 0;
2079   Int_t    qMip = 0;
2080   Int_t    nphMip = 0;
2081
2082   Float_t xTrk = 0;
2083   Float_t yTrk = 0;
2084   Float_t thetaTrk = 0;
2085   Float_t phiTrk = 0;
2086
2087   Double_t hmpPid[5]={0};
2088   Double_t hmpMom[3]={0};
2089
2090   TClonesArray &hmpidRings = *(AODEvent()->GetHMPIDrings());
2091
2092  for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
2093   {
2094     if(! esd.GetTrack(iTrack) ) continue;
2095
2096     if(esd.GetTrack(iTrack)->GetHMPIDsignal() > -20 ) {                  //
2097
2098       (esd.GetTrack(iTrack))->GetHMPIDmip(xMip, yMip, qMip, nphMip);     // Get MIP properties
2099       (esd.GetTrack(iTrack))->GetHMPIDtrk(xTrk,yTrk,thetaTrk,phiTrk);
2100       (esd.GetTrack(iTrack))->GetHMPIDpid(hmpPid);
2101       if((esd.GetTrack(iTrack))->GetOuterHmpParam()) (esd.GetTrack(iTrack))->GetOuterHmpPxPyPz(hmpMom);
2102
2103      if(esd.GetTrack(iTrack)->GetHMPIDsignal() == 0 && thetaTrk == 0 && qMip == 0 && nphMip ==0 ) continue;  //
2104
2105      new(hmpidRings[cntHmpidGoodTracks++]) AliAODHMPIDrings(
2106                                                                  (esd.GetTrack(iTrack))->GetID(),             // Unique track id to attach the ring to
2107                                                                  1000000*nphMip+qMip,                         // MIP charge and number of photons
2108                                                                  (esd.GetTrack(iTrack))->GetHMPIDcluIdx(),    // 1000000*chamber id + cluster idx of the assigned MIP cluster
2109                                                                  thetaTrk,                                    // track inclination angle theta
2110                                                                  phiTrk,                                      // track inclination angle phi
2111                                                                  (esd.GetTrack(iTrack))->GetHMPIDsignal(),    // Cherenkov angle
2112                                                                  (esd.GetTrack(iTrack))->GetHMPIDoccupancy(), // Occupancy claculated for the given chamber
2113                                                                  (esd.GetTrack(iTrack))->GetHMPIDchi2(),      // Ring resolution squared
2114                                                                   xTrk,                                       // Track x coordinate (LORS)
2115                                                                   yTrk,                                       // Track y coordinate (LORS)
2116                                                                   xMip,                                       // MIP x coordinate (LORS)
2117                                                                   yMip,                                       // MIP y coordinate (LORS)
2118                                                                   hmpPid,                                     // PID probablities from ESD, remove later once it is in CombinedPid
2119                                                                   hmpMom                                      // Track momentum in HMPID at ring reconstruction
2120                                                                );
2121
2122       //  Printf(Form("+++++++++ yes/no: %d  %lf %lf %lf %lf %lf %lf ",(esd.GetTrack(iTrack))->IsHMPID(),thetaTrk, (esd.GetTrack(iTrack))->GetHMPIDchi2(),xTrk, yTrk , xMip, yMip));
2123
2124
2125    }// HMPID signal > -20
2126   }//___esd track loop
2127
2128   return cntHmpidGoodTracks;
2129 }
2130
2131 //______________________________________________________________________________
2132 void AliAnalysisTaskESDfilter::ConvertESDtoAOD()
2133 {
2134   // ESD Filter analysis task executed for each event
2135
2136   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
2137
2138   if(!esd)return;
2139
2140   AliCodeTimerAuto("",0);
2141
2142   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );
2143
2144       // Reconstruct cascades and V0 here
2145   if (fIsV0CascadeRecoEnabled) {
2146     esd->ResetCascades();
2147     esd->ResetV0s();
2148
2149     AliV0vertexer lV0vtxer;
2150     AliCascadeVertexer lCascVtxer;
2151
2152     lV0vtxer.SetCuts(fV0Cuts);
2153     lCascVtxer.SetCuts(fCascadeCuts);
2154
2155
2156     lV0vtxer.Tracks2V0vertices(esd);
2157     lCascVtxer.V0sTracks2CascadeVertices(esd);
2158   }
2159
2160
2161   fNumberOfTracks = 0;
2162   fNumberOfPositiveTracks = 0;
2163   fNumberOfV0s = 0;
2164   fNumberOfVertices = 0;
2165   fNumberOfCascades = 0;
2166   fNumberOfKinks = 0;
2167
2168   AliAODHeader* header = ConvertHeader(*esd);
2169
2170   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);
2171   if ( fIsTZEROEnabled ) ConvertTZERO(*esd);
2172
2173   // Fetch Stack for debuggging if available
2174   fMChandler=0x0;
2175   if(MCEvent())
2176   {
2177     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
2178   }
2179
2180   // loop over events and fill them
2181   // Multiplicity information needed by the header (to be revised!)
2182   Int_t nTracks    = esd->GetNumberOfTracks();
2183   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);
2184
2185   // Update the header
2186
2187   Int_t nV0s      = esd->GetNumberOfV0s();
2188   Int_t nCascades = esd->GetNumberOfCascades();
2189   Int_t nKinks    = esd->GetNumberOfKinks();
2190   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
2191   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
2192   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
2193   nVertices+=nPileSPDVertices;
2194   nVertices+=nPileTrkVertices;
2195   Int_t nJets     = 0;
2196   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
2197   Int_t nFmdClus  = 0;
2198   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
2199   Int_t nHmpidRings = 0;
2200
2201   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));
2202
2203   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus,nHmpidRings);
2204
2205   if (nV0s > 0)
2206   {
2207     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0
2208     fAODV0VtxRefs = new TRefArray(nV0s);
2209     // RefArray to store the mapping between esd V0 number and newly created AOD-V0
2210     fAODV0Refs = new TRefArray(nV0s);
2211     // Array to take into account the V0s already added to the AOD (V0 within cascades)
2212     fUsedV0 = new Bool_t[nV0s];
2213     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;
2214   }
2215
2216   if (nTracks>0)
2217   {
2218     // RefArray to store the mapping between esd track number and newly created AOD-Track
2219
2220     fAODTrackRefs = new TRefArray(nTracks);
2221
2222     // Array to take into account the tracks already added to the AOD
2223     fUsedTrack = new Bool_t[nTracks];
2224     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;
2225   }
2226
2227   // Array to take into account the kinks already added to the AOD
2228   if (nKinks>0)
2229   {
2230     fUsedKink = new Bool_t[nKinks];
2231     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;
2232   }
2233
2234   ConvertPrimaryVertices(*esd);
2235
2236   //setting best TOF PID
2237   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
2238   if (esdH)
2239       fESDpid = esdH->GetESDpid();
2240
2241   if (fIsPidOwner && fESDpid){
2242     delete fESDpid;
2243     fESDpid = 0;
2244   }
2245   if(!fESDpid)
2246   { //in case of no Tender attached
2247     fESDpid = new AliESDpid;
2248     fIsPidOwner = kTRUE;
2249   }
2250
2251   if(!esd->GetTOFHeader())
2252   { //protection in case the pass2 LHC10b,c,d have been processed without tender.
2253     Float_t t0spread[10];
2254     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!!
2255     for (Int_t i=0; i<10; i++) t0spread[i] = (TMath::Sqrt(esd->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
2256     fESDpid->GetTOFResponse().SetT0resolution(t0spread);
2257     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
2258       //    fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);
2259     AliTOFHeader tmpTOFHeader(0,t0spread[0],0,NULL,NULL,NULL,intrinsicTOFres,t0spread[0]);
2260     AODEvent()->SetTOFHeader(&tmpTOFHeader);         // write dummy TOF header in AOD
2261   } else {
2262     AODEvent()->SetTOFHeader(esd->GetTOFHeader());    // write TOF header in AOD
2263   }
2264
2265   //  if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender.
2266
2267   if ( fAreCascadesEnabled ) ConvertCascades(*esd);
2268
2269   if ( fAreV0sEnabled ) ConvertV0s(*esd);
2270
2271   if ( fAreKinksEnabled ) ConvertKinks(*esd);
2272
2273   if ( fAreTracksEnabled ) ConvertTracks(*esd);
2274
2275   // Update number of AOD tracks in header at the end of track loop (M.G.)
2276   header->SetRefMultiplicity(fNumberOfTracks);
2277   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);
2278   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);
2279
2280   if ( fTPCConstrainedFilterMask ) ConvertTPCOnlyTracks(*esd);
2281   if( fGlobalConstrainedFilterMask) ConvertGlobalConstrainedTracks(*esd);
2282
2283   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);
2284
2285   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);
2286
2287   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);
2288
2289   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);
2290
2291         if ( fAreEMCALTriggerEnabled )ConvertCaloTrigger(TString("EMCAL"), *esd);
2292
2293         if ( fArePHOSTriggerEnabled )ConvertCaloTrigger(TString("PHOS"), *esd);
2294
2295   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);
2296   if ( fIsZDCEnabled ) ConvertZDC(*esd);
2297
2298   if(fIsHMPIDEnabled) nHmpidRings = ConvertHMPID(*esd);
2299
2300   delete fAODTrackRefs; fAODTrackRefs=0x0;
2301   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;
2302   delete fAODV0Refs; fAODV0Refs=0x0;
2303
2304   delete[] fUsedTrack; fUsedTrack=0x0;
2305   delete[] fUsedV0; fUsedV0=0x0;
2306   delete[] fUsedKink; fUsedKink=0x0;
2307
2308   if ( fIsPidOwner){
2309     delete fESDpid;
2310     fESDpid = 0x0;
2311   }
2312
2313
2314 }
2315
2316
2317 //______________________________________________________________________________
2318 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)
2319 {
2320   //
2321   // Setter for the raw PID detector signals
2322   //
2323
2324   // Save PID object for candidate electrons
2325     Bool_t pidSave = kFALSE;
2326     if (fTrackFilter) {
2327         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");
2328         if (selectInfo)  pidSave = kTRUE;
2329     }
2330
2331
2332     // Tracks passing pt cut
2333     if(esdtrack->Pt()>fHighPthreshold) {
2334         pidSave = kTRUE;
2335     } else {
2336         if(fPtshape){
2337             if(esdtrack->Pt()> fPtshape->GetXmin()){
2338                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
2339                 if(gRandom->Rndm(0)<1./y){
2340                     pidSave = kTRUE;
2341                 }//end rndm
2342             }//end if p < pmin
2343         }//end if p function
2344     }// end else
2345
2346     if (pidSave) {
2347       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track
2348         detpid = new AliAODPid();
2349         SetDetectorRawSignals(detpid,esdtrack);
2350         aodtrack->SetDetPID(detpid);
2351       }
2352     }
2353 }
2354
2355 //______________________________________________________________________________
2356 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)
2357 {
2358 //
2359 //assignment of the detector signals (AliXXXesdPID inspired)
2360 //
2361  if(!track){
2362  AliInfo("no ESD track found. .....exiting");
2363  return;
2364  }
2365  // TPC momentum
2366  const AliExternalTrackParam *in=track->GetInnerParam();
2367  if (in) {
2368    aodpid->SetTPCmomentum(in->GetP());
2369  }else{
2370    aodpid->SetTPCmomentum(-1.);
2371  }
2372
2373
2374  aodpid->SetITSsignal(track->GetITSsignal());
2375  Double_t itsdedx[4]; // dE/dx samples for individual ITS layers
2376  track->GetITSdEdxSamples(itsdedx);
2377  aodpid->SetITSdEdxSamples(itsdedx);
2378
2379  aodpid->SetTPCsignal(track->GetTPCsignal());
2380  aodpid->SetTPCsignalN(track->GetTPCsignalN());
2381  if(track->GetTPCdEdxInfo()) aodpid->SetTPCdEdxInfo(track->GetTPCdEdxInfo());
2382
2383  //n TRD planes = 6
2384  Int_t nslices = track->GetNumberOfTRDslices()*6;
2385  TArrayD trdslices(nslices);
2386  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
2387    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
2388  }
2389
2390 //TRD momentum
2391  for(Int_t iPl=0;iPl<6;iPl++){
2392    Double_t trdmom=track->GetTRDmomentum(iPl);
2393    aodpid->SetTRDmomentum(iPl,trdmom);
2394  }
2395
2396  aodpid->SetTRDsignal(track->GetNumberOfTRDslices()*6,trdslices.GetArray());
2397
2398  //TRD clusters and tracklets
2399  aodpid->SetTRDncls(track->GetTRDncls());
2400  aodpid->SetTRDntrackletsPID(track->GetTRDntrackletsPID());
2401
2402  aodpid->SetTRDChi2(track->GetTRDchi2());
2403
2404  //TOF PID
2405  Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
2406  aodpid->SetIntegratedTimes(times);
2407
2408    //  Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
2409    //  aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);
2410    aodpid->SetTOFsignal(track->GetTOFsignal());
2411
2412   Double_t tofRes[5];
2413   for (Int_t iMass=0; iMass<5; iMass++){
2414     //    tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));
2415     tofRes[iMass]=0; //backward compatibility
2416   }
2417   aodpid->SetTOFpidResolution(tofRes);
2418
2419 //  aodpid->SetHMPIDsignal(0); // set to zero for compression but it will be removed later
2420
2421 }
2422
2423 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)
2424 {
2425     // Calculate chi2 per ndf for track
2426     Int_t  nClustersTPC = track->GetTPCNcls();
2427
2428     if ( nClustersTPC > 5) {
2429        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));
2430     } else {
2431        return (-1.);
2432     }
2433  }
2434
2435
2436 //______________________________________________________________________________
2437 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
2438 {
2439 // Terminate analysis
2440 //
2441     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
2442 }
2443
2444 //______________________________________________________________________________
2445 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){
2446 // Print MC info
2447   if(!pStack)return;
2448   label = TMath::Abs(label);
2449   TParticle *part = pStack->Particle(label);
2450   Printf("########################");
2451   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
2452   part->Print();
2453   TParticle* mother = part;
2454   Int_t imo = part->GetFirstMother();
2455   Int_t nprim = pStack->GetNprimary();
2456   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
2457   while((imo >= nprim)) {
2458     mother =  pStack->Particle(imo);
2459     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
2460     mother->Print();
2461     imo =  mother->GetFirstMother();
2462   }
2463   Printf("########################");
2464 }
2465
2466 //______________________________________________________
2467