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