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