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