PWG/FLOW/Tasks/AliFlowTrackCuts.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 // AliFlowTrackCuts:
  19 // Data selection for flow framework
  20 //
  21 // origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
  22 //
  23 // This class gurantees consistency of cut methods, trackparameter
  24 // selection (global tracks, TPC only, etc..) and parameter mixing
  25 // in the flow framework. Transparently handles different input types:
  26 // ESD, MC, AOD.
  27 // This class works in 2 steps: first the requested track parameters are
  28 // constructed (to be set by SetParamType() ), then cuts are applied.
  29 // the constructed track can be requested AFTER checking the cuts by
  30 // calling GetTrack(), in this case the cut object stays in control,
  31 // caller does not have to delete the track.
  32 // Additionally caller can request an AliFlowTrack object to be constructed
  33 // according the parameter mixing scenario requested by SetParamMix().
  34 // AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
  35 // so caller needs to take care of the freshly created object.
  36
  37 #include <limits.h>
  38 #include <float.h>
  39 #include <TMatrix.h>
  40 #include "TMCProcess.h"
  41 #include "TParticle.h"
  42 #include "TH2F.h"
  43 #include "AliStack.h"
  44 #include "TBrowser.h"
  45 #include "AliMCEvent.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliVParticle.h"
  49 #include "AliMCParticle.h"
  50 #include "AliESDtrack.h"
  51 #include "AliMultiplicity.h"
  52 #include "AliAODTrack.h"
  53 #include "AliFlowTrackSimple.h"
  54 #include "AliFlowTrack.h"
  55 #include "AliFlowTrackCuts.h"
  56 #include "AliLog.h"
  57 #include "AliESDpid.h"
  58 #include "AliESDPmdTrack.h"
  59 #include "AliESDUtils.h" //TODO
  60 #include "AliFlowBayesianPID.h"
  61
  62 ClassImp(AliFlowTrackCuts)
  63
  64 //-----------------------------------------------------------------------
  65 AliFlowTrackCuts::AliFlowTrackCuts():
  66   AliFlowTrackSimpleCuts(),
  67   fAliESDtrackCuts(NULL),
  68   fQA(NULL),
  69   fCutMC(kFALSE),
  70   fCutMChasTrackReferences(kFALSE),
  71   fCutMCprocessType(kFALSE),
  72   fMCprocessType(kPNoProcess),
  73   fCutMCPID(kFALSE),
  74   fMCPID(0),
  75   fCutMCfirstMotherPID(kFALSE),
  76   fMCfirstMotherPID(0),
  77   fIgnoreSignInMCPID(kFALSE),
  78   fCutMCisPrimary(kFALSE),
  79   fRequireTransportBitForPrimaries(kTRUE),
  80   fMCisPrimary(kFALSE),
  81   fRequireCharge(kFALSE),
  82   fFakesAreOK(kTRUE),
  83   fCutSPDtrackletDeltaPhi(kFALSE),
  84   fSPDtrackletDeltaPhiMax(FLT_MAX),
  85   fSPDtrackletDeltaPhiMin(-FLT_MAX),
  86   fIgnoreTPCzRange(kFALSE),
  87   fIgnoreTPCzRangeMax(FLT_MAX),
  88   fIgnoreTPCzRangeMin(-FLT_MAX),
  89   fCutChi2PerClusterTPC(kFALSE),
  90   fMaxChi2PerClusterTPC(FLT_MAX),
  91   fMinChi2PerClusterTPC(-FLT_MAX),
  92   fCutNClustersTPC(kFALSE),
  93   fNClustersTPCMax(INT_MAX),
  94   fNClustersTPCMin(INT_MIN),
  95   fCutNClustersITS(kFALSE),
  96   fNClustersITSMax(INT_MAX),
  97   fNClustersITSMin(INT_MIN),
  98   fUseAODFilterBit(kFALSE),
  99   fAODFilterBit(0),
 100   fCutDCAToVertexXY(kFALSE),
 101   fCutDCAToVertexZ(kFALSE),
 102   fCutMinimalTPCdedx(kFALSE),
 103   fMinimalTPCdedx(0.),
 104   fLinearizeVZEROresponse(kFALSE),
 105   fCutPmdDet(kFALSE),
 106   fPmdDet(0),
 107   fCutPmdAdc(kFALSE),
 108   fPmdAdc(0.),
 109   fCutPmdNcell(kFALSE),
 110   fPmdNcell(0.),
 111   fParamType(kGlobal),
 112   fParamMix(kPure),
 113   fTrack(NULL),
 114   fTrackPhi(0.),
 115   fTrackEta(0.),
 116   fTrackWeight(1.),
 117   fTrackLabel(INT_MIN),
 118   fMCevent(NULL),
 119   fMCparticle(NULL),
 120   fEvent(NULL),
 121   fTPCtrack(),
 122   fFlowTagType(AliFlowTrackSimple::kInvalid),
 123   fESDpid(),
 124   fBayesianResponse(NULL),
 125   fPIDsource(kTOFpid),
 126   fTPCpidCuts(NULL),
 127   fTOFpidCuts(NULL),
 128   fParticleID(AliPID::kUnknown),
 129   fParticleProbability(.9),
 130   fAllowTOFmismatchFlag(kFALSE),
 131   fRequireStrictTOFTPCagreement(kFALSE),
 132   fCutRejectElectronsWithTPCpid(kFALSE),
 133   fProbBayes(0.0),
 134   fCurrCentr(0.0)
 135 {
 136   //io constructor
 137   SetPriors(); //init arrays
 138 }
 139
 140 //-----------------------------------------------------------------------
 141 AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
 142   AliFlowTrackSimpleCuts(name),
 143   fAliESDtrackCuts(NULL),
 144   fQA(NULL),
 145   fCutMC(kFALSE),
 146   fCutMChasTrackReferences(kFALSE),
 147   fCutMCprocessType(kFALSE),
 148   fMCprocessType(kPNoProcess),
 149   fCutMCPID(kFALSE),
 150   fMCPID(0),
 151   fCutMCfirstMotherPID(kFALSE),
 152   fMCfirstMotherPID(0),
 153   fIgnoreSignInMCPID(kFALSE),
 154   fCutMCisPrimary(kFALSE),
 155   fRequireTransportBitForPrimaries(kTRUE),
 156   fMCisPrimary(kFALSE),
 157   fRequireCharge(kFALSE),
 158   fFakesAreOK(kTRUE),
 159   fCutSPDtrackletDeltaPhi(kFALSE),
 160   fSPDtrackletDeltaPhiMax(FLT_MAX),
 161   fSPDtrackletDeltaPhiMin(-FLT_MAX),
 162   fIgnoreTPCzRange(kFALSE),
 163   fIgnoreTPCzRangeMax(FLT_MAX),
 164   fIgnoreTPCzRangeMin(-FLT_MAX),
 165   fCutChi2PerClusterTPC(kFALSE),
 166   fMaxChi2PerClusterTPC(FLT_MAX),
 167   fMinChi2PerClusterTPC(-FLT_MAX),
 168   fCutNClustersTPC(kFALSE),
 169   fNClustersTPCMax(INT_MAX),
 170   fNClustersTPCMin(INT_MIN),
 171   fCutNClustersITS(kFALSE),
 172   fNClustersITSMax(INT_MAX),
 173   fNClustersITSMin(INT_MIN),
 174   fUseAODFilterBit(kFALSE),
 175   fAODFilterBit(0),
 176   fCutDCAToVertexXY(kFALSE),
 177   fCutDCAToVertexZ(kFALSE),
 178   fCutMinimalTPCdedx(kFALSE),
 179   fMinimalTPCdedx(0.),
 180   fLinearizeVZEROresponse(kFALSE),
 181   fCutPmdDet(kFALSE),
 182   fPmdDet(0),
 183   fCutPmdAdc(kFALSE),
 184   fPmdAdc(0.),
 185   fCutPmdNcell(kFALSE),
 186   fPmdNcell(0.),
 187   fParamType(kGlobal),
 188   fParamMix(kPure),
 189   fTrack(NULL),
 190   fTrackPhi(0.),
 191   fTrackEta(0.),
 192   fTrackWeight(1.),
 193   fTrackLabel(INT_MIN),
 194   fMCevent(NULL),
 195   fMCparticle(NULL),
 196   fEvent(NULL),
 197   fTPCtrack(),
 198   fFlowTagType(AliFlowTrackSimple::kInvalid),
 199   fESDpid(),
 200   fBayesianResponse(NULL),
 201   fPIDsource(kTOFpid),
 202   fTPCpidCuts(NULL),
 203   fTOFpidCuts(NULL),
 204   fParticleID(AliPID::kUnknown),
 205   fParticleProbability(.9),
 206   fAllowTOFmismatchFlag(kFALSE),
 207   fRequireStrictTOFTPCagreement(kFALSE),
 208   fCutRejectElectronsWithTPCpid(kFALSE),
 209   fProbBayes(0.0),
 210   fCurrCentr(0.0)
 211 {
 212   //constructor
 213   SetTitle("AliFlowTrackCuts");
 214   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
 215                                                     2.63394e+01,
 216                                                     5.04114e-11,
 217                                                     2.12543e+00,
 218                                                     4.88663e+00 );
 219   SetPriors(); //init arrays
 220
 221   // New PID procedure (Bayesian Combined PID)
 222   fBayesianResponse = new AliFlowBayesianPID();
 223   fBayesianResponse->SetNewTrackParam();
 224 }
 225
 226 //-----------------------------------------------------------------------
 227 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
 228   AliFlowTrackSimpleCuts(that),
 229   fAliESDtrackCuts(NULL),
 230   fQA(NULL),
 231   fCutMC(that.fCutMC),
 232   fCutMChasTrackReferences(that.fCutMChasTrackReferences),
 233   fCutMCprocessType(that.fCutMCprocessType),
 234   fMCprocessType(that.fMCprocessType),
 235   fCutMCPID(that.fCutMCPID),
 236   fMCPID(that.fMCPID),
 237   fCutMCfirstMotherPID(that.fCutMCfirstMotherPID),
 238   fMCfirstMotherPID(that.fMCfirstMotherPID),
 239   fIgnoreSignInMCPID(that.fIgnoreSignInMCPID),
 240   fCutMCisPrimary(that.fCutMCisPrimary),
 241   fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
 242   fMCisPrimary(that.fMCisPrimary),
 243   fRequireCharge(that.fRequireCharge),
 244   fFakesAreOK(that.fFakesAreOK),
 245   fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
 246   fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
 247   fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
 248   fIgnoreTPCzRange(that.fIgnoreTPCzRange),
 249   fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
 250   fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
 251   fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
 252   fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
 253   fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
 254   fCutNClustersTPC(that.fCutNClustersTPC),
 255   fNClustersTPCMax(that.fNClustersTPCMax),
 256   fNClustersTPCMin(that.fNClustersTPCMin),
 257   fCutNClustersITS(that.fCutNClustersITS),
 258   fNClustersITSMax(that.fNClustersITSMax),
 259   fNClustersITSMin(that.fNClustersITSMin),
 260   fUseAODFilterBit(that.fUseAODFilterBit),
 261   fAODFilterBit(that.fAODFilterBit),
 262   fCutDCAToVertexXY(that.fCutDCAToVertexXY),
 263   fCutDCAToVertexZ(that.fCutDCAToVertexZ),
 264   fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
 265   fMinimalTPCdedx(that.fMinimalTPCdedx),
 266   fLinearizeVZEROresponse(that.fLinearizeVZEROresponse),
 267   fCutPmdDet(that.fCutPmdDet),
 268   fPmdDet(that.fPmdDet),
 269   fCutPmdAdc(that.fCutPmdAdc),
 270   fPmdAdc(that.fPmdAdc),
 271   fCutPmdNcell(that.fCutPmdNcell),
 272   fPmdNcell(that.fPmdNcell),
 273   fParamType(that.fParamType),
 274   fParamMix(that.fParamMix),
 275   fTrack(NULL),
 276   fTrackPhi(0.),
 277   fTrackEta(0.),
 278   fTrackWeight(1.),
 279   fTrackLabel(INT_MIN),
 280   fMCevent(NULL),
 281   fMCparticle(NULL),
 282   fEvent(NULL),
 283   fTPCtrack(),
 284   fFlowTagType(that.fFlowTagType),
 285   fESDpid(that.fESDpid),
 286   fBayesianResponse(NULL),
 287   fPIDsource(that.fPIDsource),
 288   fTPCpidCuts(NULL),
 289   fTOFpidCuts(NULL),
 290   fParticleID(that.fParticleID),
 291   fParticleProbability(that.fParticleProbability),
 292   fAllowTOFmismatchFlag(that.fAllowTOFmismatchFlag),
 293   fRequireStrictTOFTPCagreement(that.fRequireStrictTOFTPCagreement),
 294   fCutRejectElectronsWithTPCpid(that.fCutRejectElectronsWithTPCpid),
 295   fProbBayes(0.0),
 296   fCurrCentr(0.0)
 297 {
 298   //copy constructor
 299   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 300   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 301   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
 302   SetPriors(); //init arrays
 303   if (that.fQA) DefineHistograms();
 304
 305   // New PID procedure (Bayesian Combined PID)
 306   fBayesianResponse = new AliFlowBayesianPID();
 307   fBayesianResponse->SetNewTrackParam();
 308 }
 309
 310 //-----------------------------------------------------------------------
 311 AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
 312 {
 313   //assignment
 314   if (this==&that) return *this;
 315
 316   AliFlowTrackSimpleCuts::operator=(that);
 317   //the following may seem excessive but if AliESDtrackCuts properly does copy and clone
 318   //this approach is better memory-fragmentation-wise in some cases
 319   if (that.fAliESDtrackCuts && fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
 320   if (that.fAliESDtrackCuts && !fAliESDtrackCuts) fAliESDtrackCuts=new AliESDtrackCuts(*(that.fAliESDtrackCuts));
 321   if (!that.fAliESDtrackCuts) delete fAliESDtrackCuts; fAliESDtrackCuts=NULL;
 322   //these guys we don't need to copy, just reinit
 323   if (that.fQA) {fQA->Delete(); delete fQA; fQA=NULL; DefineHistograms();}
 324   fCutMC=that.fCutMC;
 325   fCutMChasTrackReferences=that.fCutMChasTrackReferences;
 326   fCutMCprocessType=that.fCutMCprocessType;
 327   fMCprocessType=that.fMCprocessType;
 328   fCutMCPID=that.fCutMCPID;
 329   fMCPID=that.fMCPID;
 330   fCutMCfirstMotherPID=that.fCutMCfirstMotherPID;
 331   fMCfirstMotherPID=that.fMCfirstMotherPID;
 332   fIgnoreSignInMCPID=that.fIgnoreSignInMCPID,
 333   fCutMCisPrimary=that.fCutMCisPrimary;
 334   fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
 335   fMCisPrimary=that.fMCisPrimary;
 336   fRequireCharge=that.fRequireCharge;
 337   fFakesAreOK=that.fFakesAreOK;
 338   fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
 339   fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
 340   fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
 341   fIgnoreTPCzRange=that.fIgnoreTPCzRange;
 342   fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
 343   fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
 344   fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
 345   fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
 346   fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
 347   fCutNClustersTPC=that.fCutNClustersTPC;
 348   fNClustersTPCMax=that.fNClustersTPCMax;
 349   fNClustersTPCMin=that.fNClustersTPCMin;
 350   fCutNClustersITS=that.fCutNClustersITS;
 351   fNClustersITSMax=that.fNClustersITSMax;
 352   fNClustersITSMin=that.fNClustersITSMin;
 353   fUseAODFilterBit=that.fUseAODFilterBit;
 354   fAODFilterBit=that.fAODFilterBit;
 355   fCutDCAToVertexXY=that.fCutDCAToVertexXY;
 356   fCutDCAToVertexZ=that.fCutDCAToVertexZ;
 357   fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
 358   fMinimalTPCdedx=that.fMinimalTPCdedx;
 359   fLinearizeVZEROresponse=that.fLinearizeVZEROresponse;
 360   fCutPmdDet=that.fCutPmdDet;
 361   fPmdDet=that.fPmdDet;
 362   fCutPmdAdc=that.fCutPmdAdc;
 363   fPmdAdc=that.fPmdAdc;
 364   fCutPmdNcell=that.fCutPmdNcell;
 365   fPmdNcell=that.fPmdNcell;
 366   fFlowTagType=that.fFlowTagType;
 367
 368   fParamType=that.fParamType;
 369   fParamMix=that.fParamMix;
 370
 371   fTrack=NULL;
 372   fTrackEta=0.;
 373   fTrackPhi=0.;
 374   fTrackWeight=1.;
 375   fTrackLabel=INT_MIN;
 376   fMCevent=NULL;
 377   fMCparticle=NULL;
 378   fEvent=NULL;
 379
 380   fESDpid = that.fESDpid;
 381   fPIDsource = that.fPIDsource;
 382
 383   delete fTPCpidCuts;
 384   delete fTOFpidCuts;
 385   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 386   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 387
 388   fParticleID=that.fParticleID;
 389   fParticleProbability=that.fParticleProbability;
 390   fAllowTOFmismatchFlag=that.fAllowTOFmismatchFlag;
 391   fRequireStrictTOFTPCagreement=that.fRequireStrictTOFTPCagreement;
 392   fCutRejectElectronsWithTPCpid=that.fCutRejectElectronsWithTPCpid;
 393   fProbBayes = that.fProbBayes;
 394   fCurrCentr = that.fCurrCentr;
 395
 396   // New PID procedure (Bayesian Combined PID)
 397   fBayesianResponse = new AliFlowBayesianPID();
 398   fBayesianResponse->SetNewTrackParam();
 399
 400   return *this;
 401 }
 402
 403 //-----------------------------------------------------------------------
 404 AliFlowTrackCuts::~AliFlowTrackCuts()
 405 {
 406   //dtor
 407   delete fAliESDtrackCuts;
 408   delete fTPCpidCuts;
 409   delete fTOFpidCuts;
 410   if (fQA) { fQA->SetOwner(); fQA->Delete(); delete fQA; }
 411 }
 412
 413 //-----------------------------------------------------------------------
 414 void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
 415 {
 416   //set the event
 417   Clear();
 418   fEvent=event;
 419   fMCevent=mcEvent;
 420
 421   //do the magic for ESD
 422   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
 423   if (fCutPID && myESD)
 424   {
 425     //TODO: maybe call it only for the TOF options?
 426     // Added by F. Noferini for TOF PID
 427     // old procedure now implemented inside fBayesianResponse
 428     //    fESDpid.MakePID(myESD,kFALSE);
 429     // new procedure
 430     fBayesianResponse->SetDetResponse(myESD, fCurrCentr,AliESDpid::kTOF_T0,kTRUE); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
 431     fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
 432     // End F. Noferini added part
 433   }
 434
 435   //TODO: AOD
 436 }
 437
 438 //-----------------------------------------------------------------------
 439 void AliFlowTrackCuts::SetCutMC( Bool_t b )
 440 {
 441   //will we be cutting on MC information?
 442   fCutMC=b;
 443
 444   //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
 445   if (fCutMC)
 446   {
 447     fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
 448                                                        1.75295e+01,
 449                                                        3.40030e-09,
 450                                                        1.96178e+00,
 451                                                        3.91720e+00);
 452   }
 453 }
 454
 455 //-----------------------------------------------------------------------
 456 Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
 457 {
 458   //check cuts
 459   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 460   if (vparticle) return PassesCuts(vparticle);
 461   AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
 462   if (flowtrack) return PassesCuts(flowtrack);
 463   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 464   if (tracklets) return PassesCuts(tracklets,id);
 465   AliESDPmdTrack* pmdtrack = dynamic_cast<AliESDPmdTrack*>(obj);
 466   if (pmdtrack) return PassesPMDcuts(pmdtrack);
 467   AliVEvent* vvzero = dynamic_cast<AliVEvent*>(obj); // should be removed; left for protection only
 468   if (vvzero) return PassesV0cuts(id);
 469   return kFALSE;  //default when passed wrong type of object
 470 }
 471
 472 //-----------------------------------------------------------------------
 473 Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
 474 {
 475   //check cuts
 476   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 477   if (vparticle)
 478   {
 479     return PassesMCcuts(fMCevent,(fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel());
 480   }
 481   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 482   if (tracklets)
 483   {
 484     Int_t label0 = tracklets->GetLabel(id,0);
 485     Int_t label1 = tracklets->GetLabel(id,1);
 486     Int_t label = (fFakesAreOK)?TMath::Abs(label0):label0;
 487     if (label0!=label1) label = -666;
 488     return PassesMCcuts(fMCevent,label);
 489   }
 490   return kFALSE;  //default when passed wrong type of object
 491 }
 492
 493 //-----------------------------------------------------------------------
 494 Bool_t AliFlowTrackCuts::PassesCuts(const AliFlowTrackSimple* track)
 495 {
 496   //check cuts on a flowtracksimple
 497
 498   //clean up from last iteration
 499   fTrack = NULL;
 500   return AliFlowTrackSimpleCuts::PassesCuts(track);
 501 }
 502
 503 //-----------------------------------------------------------------------
 504 Bool_t AliFlowTrackCuts::PassesCuts(const AliMultiplicity* tracklet, Int_t id)
 505 {
 506   //check cuts on a tracklets
 507
 508   if (id<0) return kFALSE;
 509
 510   //clean up from last iteration, and init label
 511   fTrack = NULL;
 512   fMCparticle=NULL;
 513   fTrackLabel=-999;
 514
 515   fTrackPhi = tracklet->GetPhi(id);
 516   fTrackEta = tracklet->GetEta(id);
 517   fTrackWeight = 1.0;
 518   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
 519   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
 520
 521   //check MC info if available
 522   //if the 2 clusters have different label track cannot be good
 523   //and should therefore not pass the mc cuts
 524   Int_t label0 = tracklet->GetLabel(id,0);
 525   Int_t label1 = tracklet->GetLabel(id,1);
 526   //if possible get label and mcparticle
 527   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
 528   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
 529   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 530   //check MC cuts
 531   if (fCutMC && !PassesMCcuts()) return kFALSE;
 532   return kTRUE;
 533 }
 534
 535 //-----------------------------------------------------------------------
 536 Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
 537 {
 538   //check the MC info
 539   if (!mcEvent) return kFALSE;
 540   if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 541   AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
 542   if (!mcparticle) {AliError("no MC track"); return kFALSE;}
 543
 544   if (fCutMCisPrimary)
 545   {
 546     if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
 547   }
 548   if (fCutMCPID)
 549   {
 550     Int_t pdgCode = mcparticle->PdgCode();
 551     if (fIgnoreSignInMCPID)
 552     {
 553       if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
 554     }
 555     else
 556     {
 557       if (fMCPID != pdgCode) return kFALSE;
 558     }
 559   }
 560   if (fCutMCfirstMotherPID)
 561   {
 562
 563     TParticle* tparticle=mcparticle->Particle();
 564     Int_t firstMotherLabel = 0;
 565     if (tparticle) { firstMotherLabel = tparticle->GetFirstMother(); }
 566     AliVParticle* firstMotherParticle = mcEvent->GetTrack(firstMotherLabel);
 567     Int_t pdgcodeFirstMother = 0;
 568     if (firstMotherParticle) { pdgcodeFirstMother = firstMotherParticle->PdgCode(); }
 569     if (pdgcodeFirstMother != fMCfirstMotherPID) return kFALSE;
 570   }
 571   if ( fCutMCprocessType )
 572   {
 573     TParticle* particle = mcparticle->Particle();
 574     Int_t processID = particle->GetUniqueID();
 575     if (processID != fMCprocessType ) return kFALSE;
 576   }
 577   if (fCutMChasTrackReferences)
 578   {
 579     if (mcparticle->GetNumberOfTrackReferences()<1) return kFALSE;
 580   }
 581   return kTRUE;
 582 }
 583
 584 //-----------------------------------------------------------------------
 585 Bool_t AliFlowTrackCuts::PassesMCcuts()
 586 {
 587   //check MC info
 588   if (!fMCevent) return kFALSE;
 589   if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 590   fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 591   return PassesMCcuts(fMCevent,fTrackLabel);
 592 }
 593
 594 //-----------------------------------------------------------------------
 595 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
 596 {
 597   //check cuts for an ESD vparticle
 598
 599   ////////////////////////////////////////////////////////////////
 600   //  start by preparing the track parameters to cut on //////////
 601   ////////////////////////////////////////////////////////////////
 602   //clean up from last iteration
 603   fTrack=NULL;
 604
 605   //get the label and the mc particle
 606   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
 607   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 608   else fMCparticle=NULL;
 609
 610   Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
 611   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
 612   AliAODTrack* aodTrack = NULL;
 613   if (esdTrack)
 614   {
 615     //for an ESD track we do some magic sometimes like constructing TPC only parameters
 616     //or doing some hybrid, handle that here
 617     HandleESDtrack(esdTrack);
 618   }
 619   else
 620   {
 621     HandleVParticle(vparticle);
 622     //now check if produced particle is MC
 623     isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
 624     aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
 625   }
 626   ////////////////////////////////////////////////////////////////
 627   ////////////////////////////////////////////////////////////////
 628
 629   if (!fTrack) return kFALSE;
 630   //because it may be different from global, not needed for aodTrack because we dont do anything funky there
 631   if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
 632
 633   Bool_t pass=kTRUE;
 634   //check the common cuts for the current particle fTrack (MC,AOD,ESD)
 635   Double_t pt = fTrack->Pt();
 636   Double_t p = fTrack->P();
 637   if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
 638   if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
 639   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
 640   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
 641   if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
 642   if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
 643   if (fCutCharge && isMCparticle)
 644   {
 645     //in case of an MC particle the charge is stored in units of 1/3|e|
 646     Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
 647     if (charge!=fCharge) pass=kFALSE;
 648   }
 649   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
 650
 651   //when additionally MC info is required
 652   if (fCutMC && !PassesMCcuts()) pass=kFALSE;
 653
 654   //the case of ESD or AOD
 655   if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
 656   if (aodTrack) { if (!PassesAODcuts(aodTrack)) { pass=kFALSE; } }
 657
 658   if (fQA)
 659   {
 660     if (fMCparticle)
 661     {
 662       TParticle* tparticle=fMCparticle->Particle();
 663       Int_t processID = tparticle->GetUniqueID();
 664       Int_t firstMotherLabel = tparticle->GetFirstMother();
 665       Bool_t primary = IsPhysicalPrimary();
 666       //TLorentzVector v;
 667       //mcparticle->Particle()->ProductionVertex(v);
 668       //Double_t prodvtxX = v.X();
 669       //Double_t prodvtxY = v.Y();
 670
 671       Int_t pdgcode = fMCparticle->PdgCode();
 672       AliVParticle* firstMotherParticle = fMCevent->GetTrack(firstMotherLabel);
 673       Int_t pdgcodeFirstMother = 0;
 674       if (firstMotherParticle) {pdgcodeFirstMother = firstMotherParticle->PdgCode();}
 675
 676       Float_t pdg = 0;
 677       switch (TMath::Abs(pdgcode))
 678       {
 679         case 11:
 680           pdg = AliPID::kElectron + 0.5; break;
 681         case 13:
 682           pdg = AliPID::kMuon + 0.5; break;
 683         case 211:
 684           pdg = AliPID::kPion + 0.5; break;
 685         case 321:
 686           pdg = AliPID::kKaon + 0.5; break;
 687         case 2212:
 688           pdg = AliPID::kProton + 0.5; break;
 689         default:
 690           pdg = AliPID::kUnknown + 0.5; break;
 691       }
 692       pdg = TMath::Sign(pdg,static_cast<Float_t>(pdgcode));
 693
 694       Float_t geantCode = 0;
 695       switch (pdgcodeFirstMother)
 696       {
 697         case 22: //#gamma
 698           geantCode=1;
 699           break;
 700         case -11: //e^{+}
 701           geantCode=2;
 702           break;
 703         case 11: //e^{-}
 704           geantCode=3;
 705           break;
 706         case 12: case 14: case 16: //#nu
 707           geantCode=4;
 708           break;
 709         case -13:
 710           geantCode=5; //#mu^{+}
 711           break;
 712         case 13:
 713           geantCode=6; //#mu^{-}
 714           break;
 715         case 111: //#pi^{0}
 716           geantCode=7;
 717           break;
 718         case 211: //#pi^{+}
 719           geantCode=8;
 720           break;
 721         case -211: //#pi^{-}
 722           geantCode=9;
 723           break;
 724         case 130: //K^{0}_{L}
 725           geantCode=10;
 726           break;
 727         case 321: //K^{+}
 728           geantCode=11;
 729           break;
 730         case -321: //K^{-}
 731           geantCode=12;
 732           break;
 733         case 2112:
 734           geantCode = 13; //n
 735           break;
 736         case 2212:
 737           geantCode = 14; //p
 738           break;
 739         case -2212:
 740           geantCode=15; //#bar{p}
 741           break;
 742         case 310:
 743           geantCode=16; //K^{0}_{S}
 744           break;
 745         case 221:
 746           geantCode=17; //#eta
 747           break;
 748         case 3122:
 749           geantCode=18; //#Lambda
 750           break;
 751         case 3222:
 752           geantCode=19; //#Sigma^{+}
 753           break;
 754         case 3212:
 755           geantCode=20; //#Sigma^{0}
 756           break;
 757         case 3112:
 758           geantCode=21; //#Sigma^{-}
 759           break;
 760         case 3322:
 761           geantCode=22; //#Theta^{0}
 762           break;
 763         case 3312:
 764           geantCode=23; //#Theta^{-}
 765           break;
 766         case 3332:
 767           geantCode=24; //#Omega^{-}
 768           break;
 769         case -2112:
 770           geantCode=25; //#bar{n}
 771           break;
 772         case -3122:
 773           geantCode=26; //#bar{#Lambda}
 774           break;
 775         case -3112:
 776           geantCode=27; //#bar{#Sigma}^{-}
 777           break;
 778         case -3212:
 779           geantCode=28; //#bar{#Sigma}^{0}
 780           break;
 781         case -3222:
 782           geantCode=29; //#bar{#Sigma}^{+}
 783           break;
 784         case -3322:
 785           geantCode=30; //#bar{#Theta}^{0}
 786           break;
 787         case -3312:
 788           geantCode=31; //#bar{#Theta}^{+}
 789           break;
 790         case -3332:
 791           geantCode=32; //#bar{#Omega}^{+}
 792           break;
 793         case -15:
 794           geantCode=33; //#tau^{+}
 795           break;
 796         case 15:
 797           geantCode=34; //#tau^{-}
 798           break;
 799         case 411:
 800           geantCode=35; //D^{+}
 801           break;
 802         case -411:
 803           geantCode=36; //D^{-}
 804           break;
 805         case 421:
 806           geantCode=37; //D^{0}
 807           break;
 808         case -421:
 809           geantCode=38; //#bar{D}^{0}
 810           break;
 811         case 431:
 812           geantCode=39; //D_{s}^{+}
 813           break;
 814         case -431:
 815           geantCode=40; //#bar{D_{s}}^{-}
 816           break;
 817         case 4122:
 818           geantCode=41; //#Lambda_{c}^{+}
 819           break;
 820         case 24:
 821           geantCode=42; //W^{+}
 822           break;
 823         case -24:
 824           geantCode=43; //W^{-}
 825           break;
 826         case 23:
 827           geantCode=44; //Z^{0}
 828           break;
 829         default:
 830           geantCode = 100;
 831           break;
 832       }
 833
 834       QAbefore(2)->Fill(p,pdg);
 835       QAbefore(3)->Fill(p,primary?0.5:-0.5);
 836       QAbefore(4)->Fill(p,static_cast<Float_t>(processID));
 837       QAbefore(7)->Fill(p,geantCode+0.5);
 838       if (pass) QAafter(2)->Fill(p,pdg);
 839       if (pass) QAafter(3)->Fill(p,primary?0.5:-0.5);
 840       if (pass) QAafter(4)->Fill(p,static_cast<Float_t>(processID));
 841       if (pass) QAafter(7)->Fill(p,geantCode);
 842     }
 843   }
 844
 845   //true by default, if we didn't set any cuts
 846   return pass;
 847 }
 848
 849 //_______________________________________________________________________
 850 Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track)
 851 {
 852   //check cuts for AOD
 853   Bool_t pass = kTRUE;
 854
 855   if (fCutNClustersTPC)
 856   {
 857     Int_t ntpccls = track->GetTPCNcls();
 858     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
 859   }
 860
 861   if (fCutNClustersITS)
 862   {
 863     Int_t nitscls = track->GetITSNcls();
 864     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
 865   }
 866
 867    if (fCutChi2PerClusterTPC)
 868   {
 869     Double_t chi2tpc = track->Chi2perNDF();
 870     if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
 871   }
 872
 873   if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
 874   if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
 875
 876   if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
 877
 878   if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
 879
 880   if (fCutDCAToVertexZ && track->ZAtDCA()>GetMaxDCAToVertexZ()) pass=kFALSE;
 881
 882   if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
 883
 884   return pass;
 885 }
 886
 887 //_______________________________________________________________________
 888 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
 889 {
 890   //check cuts on ESD tracks
 891   Bool_t pass=kTRUE;
 892   Float_t dcaxy=0.0;
 893   Float_t dcaz=0.0;
 894   track->GetImpactParameters(dcaxy,dcaz);
 895   const AliExternalTrackParam* pout = track->GetOuterParam();
 896   const AliExternalTrackParam* pin = track->GetInnerParam();
 897   if (fIgnoreTPCzRange)
 898   {
 899     if (pin&&pout)
 900     {
 901       Double_t zin = pin->GetZ();
 902       Double_t zout = pout->GetZ();
 903       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
 904       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
 905       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
 906     }
 907   }
 908
 909   Int_t ntpccls = ( fParamType==kTPCstandalone )?
 910                     track->GetTPCNclsIter1():track->GetTPCNcls();
 911   if (fCutChi2PerClusterTPC)
 912   {
 913     Float_t tpcchi2 = (fParamType==kTPCstandalone)?
 914                        track->GetTPCchi2Iter1():track->GetTPCchi2();
 915     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
 916     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
 917       pass=kFALSE;
 918   }
 919
 920   if (fCutMinimalTPCdedx)
 921   {
 922     if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
 923   }
 924
 925   if (fCutNClustersTPC)
 926   {
 927     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
 928   }
 929
 930   Int_t nitscls = track->GetNcls(0);
 931   if (fCutNClustersITS)
 932   {
 933     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
 934   }
 935
 936   //some stuff is still handled by AliESDtrackCuts class - delegate
 937   if (fAliESDtrackCuts)
 938   {
 939     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
 940   }
 941
 942   //PID part with pid QA
 943   Double_t beta = GetBeta(track);
 944   Double_t dedx = Getdedx(track);
 945   if (fQA)
 946   {
 947     if (pass) QAbefore(0)->Fill(track->GetP(),beta);
 948     if (pass) QAbefore(1)->Fill(pin->GetP(),dedx);
 949   }
 950   if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
 951   {
 952     if (!PassesESDpidCut(track)) pass=kFALSE;
 953   }
 954   if (fCutRejectElectronsWithTPCpid)
 955   {
 956     //reject electrons using standard TPC pid
 957     //TODO this should be rethought....
 958     Double_t pidTPC[AliPID::kSPECIES];
 959     track->GetTPCpid(pidTPC);
 960     if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
 961   }
 962   if (fQA)
 963   {
 964     if (pass) QAafter(0)->Fill(track->GetP(),beta);
 965     if (pass) QAafter(1)->Fill(pin->GetP(),dedx);
 966   }
 967   //end pid part with qa
 968
 969   if (fQA)
 970   {
 971     Double_t pt = track->Pt();
 972     QAbefore(5)->Fill(pt,dcaxy);
 973     QAbefore(6)->Fill(pt,dcaz);
 974     if (pass) QAafter(5)->Fill(pt,dcaxy);
 975     if (pass) QAafter(6)->Fill(pt,dcaz);
 976   }
 977
 978   return pass;
 979 }
 980
 981 //-----------------------------------------------------------------------
 982 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
 983 {
 984   //handle the general case
 985   switch (fParamType)
 986   {
 987     default:
 988       fTrack = track;
 989       break;
 990   }
 991 }
 992
 993 //-----------------------------------------------------------------------
 994 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
 995 {
 996   //handle esd track
 997   switch (fParamType)
 998   {
 999     case kGlobal:
1000       fTrack = track;
1001       break;
1002     case kTPCstandalone:
1003       if (!track->FillTPCOnlyTrack(fTPCtrack))
1004       {
1005         fTrack=NULL;
1006         fMCparticle=NULL;
1007         fTrackLabel=-998;
1008         return;
1009       }
1010       fTrack = &fTPCtrack;
1011       //recalculate the label and mc particle, they may differ as TPClabel != global label
1012       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1013       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1014       else fMCparticle=NULL;
1015       break;
1016     default:
1017       fTrack = track;
1018       break;
1019   }
1020 }
1021
1022 //-----------------------------------------------------------------------
1023 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
1024 {
1025   //calculate the number of track in given event.
1026   //if argument is NULL(default) take the event attached
1027   //by SetEvent()
1028   Int_t multiplicity = 0;
1029   if (!event)
1030   {
1031     for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
1032     {
1033       if (IsSelected(GetInputObject(i))) multiplicity++;
1034     }
1035   }
1036   else
1037   {
1038     for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
1039     {
1040       if (IsSelected(event->GetTrack(i))) multiplicity++;
1041     }
1042   }
1043   return multiplicity;
1044 }
1045
1046 //-----------------------------------------------------------------------
1047 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts()
1048 {
1049   //get standard V0 cuts
1050   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts");
1051   cuts->SetParamType(kV0);
1052   cuts->SetEtaRange( -10, +10 );
1053   cuts->SetPhiMin( 0 );
1054   cuts->SetPhiMax( TMath::TwoPi() );
1055   return cuts;
1056 }
1057
1058 //-----------------------------------------------------------------------
1059 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
1060 {
1061   //get standard cuts
1062   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global tracks");
1063   cuts->SetParamType(kGlobal);
1064   cuts->SetPtRange(0.2,5.);
1065   cuts->SetEtaRange(-0.8,0.8);
1066   cuts->SetMinNClustersTPC(70);
1067   cuts->SetMinChi2PerClusterTPC(0.1);
1068   cuts->SetMaxChi2PerClusterTPC(4.0);
1069   cuts->SetMinNClustersITS(2);
1070   cuts->SetRequireITSRefit(kTRUE);
1071   cuts->SetRequireTPCRefit(kTRUE);
1072   cuts->SetMaxDCAToVertexXY(0.3);
1073   cuts->SetMaxDCAToVertexZ(0.3);
1074   cuts->SetAcceptKinkDaughters(kFALSE);
1075   cuts->SetMinimalTPCdedx(10.);
1076
1077   return cuts;
1078 }
1079
1080 //-----------------------------------------------------------------------
1081 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()
1082 {
1083   //get standard cuts
1084   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone 2010");
1085   cuts->SetParamType(kTPCstandalone);
1086   cuts->SetPtRange(0.2,5.);
1087   cuts->SetEtaRange(-0.8,0.8);
1088   cuts->SetMinNClustersTPC(70);
1089   cuts->SetMinChi2PerClusterTPC(0.2);
1090   cuts->SetMaxChi2PerClusterTPC(4.0);
1091   cuts->SetMaxDCAToVertexXY(3.0);
1092   cuts->SetMaxDCAToVertexZ(3.0);
1093   cuts->SetDCAToVertex2D(kTRUE);
1094   cuts->SetAcceptKinkDaughters(kFALSE);
1095   cuts->SetMinimalTPCdedx(10.);
1096
1097   return cuts;
1098 }
1099
1100 //-----------------------------------------------------------------------
1101 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts()
1102 {
1103   //get standard cuts
1104   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone");
1105   cuts->SetParamType(kTPCstandalone);
1106   cuts->SetPtRange(0.2,5.);
1107   cuts->SetEtaRange(-0.8,0.8);
1108   cuts->SetMinNClustersTPC(70);
1109   cuts->SetMinChi2PerClusterTPC(0.2);
1110   cuts->SetMaxChi2PerClusterTPC(4.0);
1111   cuts->SetMaxDCAToVertexXY(3.0);
1112   cuts->SetMaxDCAToVertexZ(3.0);
1113   cuts->SetDCAToVertex2D(kTRUE);
1114   cuts->SetAcceptKinkDaughters(kFALSE);
1115   cuts->SetMinimalTPCdedx(10.);
1116
1117   return cuts;
1118 }
1119
1120 //-----------------------------------------------------------------------
1121 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
1122 {
1123   //get standard cuts
1124   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
1125   delete cuts->fAliESDtrackCuts;
1126   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
1127   cuts->SetParamType(kGlobal);
1128   return cuts;
1129 }
1130
1131 //-----------------------------------------------------------------------
1132 Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
1133 {
1134   //fill a flow track from tracklet,vzero,pmd,...
1135   TParticle *tmpTParticle=NULL;
1136   AliMCParticle* tmpAliMCParticle=NULL;
1137   switch (fParamMix)
1138   {
1139     case kPure:
1140       flowtrack->SetPhi(fTrackPhi);
1141       flowtrack->SetEta(fTrackEta);
1142       flowtrack->SetWeight(fTrackWeight);
1143       break;
1144     case kTrackWithMCkine:
1145       if (!fMCparticle) return kFALSE;
1146       flowtrack->SetPhi( fMCparticle->Phi() );
1147       flowtrack->SetEta( fMCparticle->Eta() );
1148       flowtrack->SetPt( fMCparticle->Pt() );
1149       break;
1150     case kTrackWithMCpt:
1151       if (!fMCparticle) return kFALSE;
1152       flowtrack->SetPhi(fTrackPhi);
1153       flowtrack->SetEta(fTrackEta);
1154       flowtrack->SetWeight(fTrackWeight);
1155       flowtrack->SetPt(fMCparticle->Pt());
1156       break;
1157     case kTrackWithPtFromFirstMother:
1158       if (!fMCparticle) return kFALSE;
1159       flowtrack->SetPhi(fTrackPhi);
1160       flowtrack->SetEta(fTrackEta);
1161       flowtrack->SetWeight(fTrackWeight);
1162       tmpTParticle = fMCparticle->Particle();
1163       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1164       flowtrack->SetPt(tmpAliMCParticle->Pt());
1165       break;
1166     default:
1167       flowtrack->SetPhi(fTrackPhi);
1168       flowtrack->SetEta(fTrackEta);
1169       flowtrack->SetWeight(fTrackWeight);
1170       break;
1171   }
1172   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1173   return kTRUE;
1174 }
1175
1176 //-----------------------------------------------------------------------
1177 Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
1178 {
1179   //fill flow track from AliVParticle (ESD,AOD,MC)
1180   if (!fTrack) return kFALSE;
1181   TParticle *tmpTParticle=NULL;
1182   AliMCParticle* tmpAliMCParticle=NULL;
1183   AliExternalTrackParam* externalParams=NULL;
1184   AliESDtrack* esdtrack=NULL;
1185   switch(fParamMix)
1186   {
1187     case kPure:
1188       flowtrack->Set(fTrack);
1189       break;
1190     case kTrackWithMCkine:
1191       flowtrack->Set(fMCparticle);
1192       break;
1193     case kTrackWithMCPID:
1194       flowtrack->Set(fTrack);
1195       //flowtrack->setPID(...) from mc, when implemented
1196       break;
1197     case kTrackWithMCpt:
1198       if (!fMCparticle) return kFALSE;
1199       flowtrack->Set(fTrack);
1200       flowtrack->SetPt(fMCparticle->Pt());
1201       break;
1202     case kTrackWithPtFromFirstMother:
1203       if (!fMCparticle) return kFALSE;
1204       flowtrack->Set(fTrack);
1205       tmpTParticle = fMCparticle->Particle();
1206       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1207       flowtrack->SetPt(tmpAliMCParticle->Pt());
1208       break;
1209     case kTrackWithTPCInnerParams:
1210       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1211       if (!esdtrack) return kFALSE;
1212       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1213       if (!externalParams) return kFALSE;
1214       flowtrack->Set(externalParams);
1215       break;
1216     default:
1217       flowtrack->Set(fTrack);
1218       break;
1219   }
1220   if (fParamType==kMC)
1221   {
1222     flowtrack->SetSource(AliFlowTrack::kFromMC);
1223     flowtrack->SetID(fTrackLabel);
1224   }
1225   else if (dynamic_cast<AliESDtrack*>(fTrack))
1226   {
1227     flowtrack->SetSource(AliFlowTrack::kFromESD);
1228     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1229   }
1230   else if (dynamic_cast<AliAODTrack*>(fTrack))
1231   {
1232     flowtrack->SetSource(AliFlowTrack::kFromAOD);
1233     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1234   }
1235   else if (dynamic_cast<AliMCParticle*>(fTrack))
1236   {
1237     flowtrack->SetSource(AliFlowTrack::kFromMC);
1238     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1239   }
1240   return kTRUE;
1241 }
1242
1243 //-----------------------------------------------------------------------
1244 Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
1245 {
1246   //fill a flow track constructed from whatever we applied cuts on
1247   //return true on success
1248   switch (fParamType)
1249   {
1250     case kSPDtracklet:
1251       return FillFlowTrackGeneric(track);
1252     case kPMD:
1253       return FillFlowTrackGeneric(track);
1254     case kV0:
1255       return FillFlowTrackGeneric(track);
1256     default:
1257       return FillFlowTrackVParticle(track);
1258   }
1259 }
1260
1261 //-----------------------------------------------------------------------
1262 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
1263 {
1264   //make a flow track from tracklet
1265   AliFlowTrack* flowtrack=NULL;
1266   TParticle *tmpTParticle=NULL;
1267   AliMCParticle* tmpAliMCParticle=NULL;
1268   switch (fParamMix)
1269   {
1270     case kPure:
1271       flowtrack = new AliFlowTrack();
1272       flowtrack->SetPhi(fTrackPhi);
1273       flowtrack->SetEta(fTrackEta);
1274       flowtrack->SetWeight(fTrackWeight);
1275       break;
1276     case kTrackWithMCkine:
1277       if (!fMCparticle) return NULL;
1278       flowtrack = new AliFlowTrack();
1279       flowtrack->SetPhi( fMCparticle->Phi() );
1280       flowtrack->SetEta( fMCparticle->Eta() );
1281       flowtrack->SetPt( fMCparticle->Pt() );
1282       break;
1283     case kTrackWithMCpt:
1284       if (!fMCparticle) return NULL;
1285       flowtrack = new AliFlowTrack();
1286       flowtrack->SetPhi(fTrackPhi);
1287       flowtrack->SetEta(fTrackEta);
1288       flowtrack->SetWeight(fTrackWeight);
1289       flowtrack->SetPt(fMCparticle->Pt());
1290       break;
1291     case kTrackWithPtFromFirstMother:
1292       if (!fMCparticle) return NULL;
1293       flowtrack = new AliFlowTrack();
1294       flowtrack->SetPhi(fTrackPhi);
1295       flowtrack->SetEta(fTrackEta);
1296       flowtrack->SetWeight(fTrackWeight);
1297       tmpTParticle = fMCparticle->Particle();
1298       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1299       flowtrack->SetPt(tmpAliMCParticle->Pt());
1300       break;
1301     default:
1302       flowtrack = new AliFlowTrack();
1303       flowtrack->SetPhi(fTrackPhi);
1304       flowtrack->SetEta(fTrackEta);
1305       flowtrack->SetWeight(fTrackWeight);
1306       break;
1307   }
1308   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1309   return flowtrack;
1310 }
1311
1312 //-----------------------------------------------------------------------
1313 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
1314 {
1315   //make flow track from AliVParticle (ESD,AOD,MC)
1316   if (!fTrack) return NULL;
1317   AliFlowTrack* flowtrack=NULL;
1318   TParticle *tmpTParticle=NULL;
1319   AliMCParticle* tmpAliMCParticle=NULL;
1320   AliExternalTrackParam* externalParams=NULL;
1321   AliESDtrack* esdtrack=NULL;
1322   switch(fParamMix)
1323   {
1324     case kPure:
1325       flowtrack = new AliFlowTrack(fTrack);
1326       break;
1327     case kTrackWithMCkine:
1328       flowtrack = new AliFlowTrack(fMCparticle);
1329       break;
1330     case kTrackWithMCPID:
1331       flowtrack = new AliFlowTrack(fTrack);
1332       //flowtrack->setPID(...) from mc, when implemented
1333       break;
1334     case kTrackWithMCpt:
1335       if (!fMCparticle) return NULL;
1336       flowtrack = new AliFlowTrack(fTrack);
1337       flowtrack->SetPt(fMCparticle->Pt());
1338       break;
1339     case kTrackWithPtFromFirstMother:
1340       if (!fMCparticle) return NULL;
1341       flowtrack = new AliFlowTrack(fTrack);
1342       tmpTParticle = fMCparticle->Particle();
1343       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1344       flowtrack->SetPt(tmpAliMCParticle->Pt());
1345       break;
1346     case kTrackWithTPCInnerParams:
1347       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1348       if (!esdtrack) return NULL;
1349       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1350       if (!externalParams) return NULL;
1351       flowtrack = new AliFlowTrack(externalParams);
1352       break;
1353     default:
1354       flowtrack = new AliFlowTrack(fTrack);
1355       break;
1356   }
1357   if (fParamType==kMC)
1358   {
1359     flowtrack->SetSource(AliFlowTrack::kFromMC);
1360     flowtrack->SetID(fTrackLabel);
1361   }
1362   else if (dynamic_cast<AliESDtrack*>(fTrack))
1363   {
1364     flowtrack->SetSource(AliFlowTrack::kFromESD);
1365     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1366   }
1367   else if (dynamic_cast<AliAODTrack*>(fTrack))
1368   {
1369     flowtrack->SetSource(AliFlowTrack::kFromAOD);
1370     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1371   }
1372   else if (dynamic_cast<AliMCParticle*>(fTrack))
1373   {
1374     flowtrack->SetSource(AliFlowTrack::kFromMC);
1375     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1376   }
1377   return flowtrack;
1378 }
1379
1380 //-----------------------------------------------------------------------
1381 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
1382 {
1383   //make a flow track from PMD track
1384   AliFlowTrack* flowtrack=NULL;
1385   TParticle *tmpTParticle=NULL;
1386   AliMCParticle* tmpAliMCParticle=NULL;
1387   switch (fParamMix)
1388   {
1389     case kPure:
1390       flowtrack = new AliFlowTrack();
1391       flowtrack->SetPhi(fTrackPhi);
1392       flowtrack->SetEta(fTrackEta);
1393       flowtrack->SetWeight(fTrackWeight);
1394       break;
1395     case kTrackWithMCkine:
1396       if (!fMCparticle) return NULL;
1397       flowtrack = new AliFlowTrack();
1398       flowtrack->SetPhi( fMCparticle->Phi() );
1399       flowtrack->SetEta( fMCparticle->Eta() );
1400       flowtrack->SetWeight(fTrackWeight);
1401       flowtrack->SetPt( fMCparticle->Pt() );
1402       break;
1403     case kTrackWithMCpt:
1404       if (!fMCparticle) return NULL;
1405       flowtrack = new AliFlowTrack();
1406       flowtrack->SetPhi(fTrackPhi);
1407       flowtrack->SetEta(fTrackEta);
1408       flowtrack->SetWeight(fTrackWeight);
1409       flowtrack->SetPt(fMCparticle->Pt());
1410       break;
1411     case kTrackWithPtFromFirstMother:
1412       if (!fMCparticle) return NULL;
1413       flowtrack = new AliFlowTrack();
1414       flowtrack->SetPhi(fTrackPhi);
1415       flowtrack->SetEta(fTrackEta);
1416       flowtrack->SetWeight(fTrackWeight);
1417       tmpTParticle = fMCparticle->Particle();
1418       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1419       flowtrack->SetPt(tmpAliMCParticle->Pt());
1420       break;
1421     default:
1422       flowtrack = new AliFlowTrack();
1423       flowtrack->SetPhi(fTrackPhi);
1424       flowtrack->SetEta(fTrackEta);
1425       flowtrack->SetWeight(fTrackWeight);
1426       break;
1427   }
1428
1429   flowtrack->SetSource(AliFlowTrack::kFromPMD);
1430   return flowtrack;
1431 }
1432
1433 //-----------------------------------------------------------------------
1434 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackV0() const
1435 {
1436   //make a flow track from V0
1437   AliFlowTrack* flowtrack=NULL;
1438   TParticle *tmpTParticle=NULL;
1439   AliMCParticle* tmpAliMCParticle=NULL;
1440   switch (fParamMix)
1441   {
1442     case kPure:
1443       flowtrack = new AliFlowTrack();
1444       flowtrack->SetPhi(fTrackPhi);
1445       flowtrack->SetEta(fTrackEta);
1446       flowtrack->SetWeight(fTrackWeight);
1447       break;
1448     case kTrackWithMCkine:
1449       if (!fMCparticle) return NULL;
1450       flowtrack = new AliFlowTrack();
1451       flowtrack->SetPhi( fMCparticle->Phi() );
1452       flowtrack->SetEta( fMCparticle->Eta() );
1453       flowtrack->SetWeight(fTrackWeight);
1454       flowtrack->SetPt( fMCparticle->Pt() );
1455       break;
1456     case kTrackWithMCpt:
1457       if (!fMCparticle) return NULL;
1458       flowtrack = new AliFlowTrack();
1459       flowtrack->SetPhi(fTrackPhi);
1460       flowtrack->SetEta(fTrackEta);
1461       flowtrack->SetWeight(fTrackWeight);
1462       flowtrack->SetPt(fMCparticle->Pt());
1463       break;
1464     case kTrackWithPtFromFirstMother:
1465       if (!fMCparticle) return NULL;
1466       flowtrack = new AliFlowTrack();
1467       flowtrack->SetPhi(fTrackPhi);
1468       flowtrack->SetEta(fTrackEta);
1469       flowtrack->SetWeight(fTrackWeight);
1470       tmpTParticle = fMCparticle->Particle();
1471       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1472       flowtrack->SetPt(tmpAliMCParticle->Pt());
1473       break;
1474     default:
1475       flowtrack = new AliFlowTrack();
1476       flowtrack->SetPhi(fTrackPhi);
1477       flowtrack->SetEta(fTrackEta);
1478       flowtrack->SetWeight(fTrackWeight);
1479       break;
1480   }
1481
1482   flowtrack->SetSource(AliFlowTrack::kFromV0);
1483   return flowtrack;
1484 }
1485
1486 //-----------------------------------------------------------------------
1487 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
1488 {
1489   //get a flow track constructed from whatever we applied cuts on
1490   //caller is resposible for deletion
1491   //if construction fails return NULL
1492   //TODO: for tracklets, PMD and V0 we probably need just one method,
1493   //something like MakeFlowTrackGeneric(), wait with this until
1494   //requirements quirks are known.
1495   switch (fParamType)
1496   {
1497     case kSPDtracklet:
1498       return MakeFlowTrackSPDtracklet();
1499     case kPMD:
1500       return MakeFlowTrackPMDtrack();
1501     case kV0:
1502       return MakeFlowTrackV0();
1503     default:
1504       return MakeFlowTrackVParticle();
1505   }
1506 }
1507
1508 //-----------------------------------------------------------------------
1509 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
1510 {
1511   //check if current particle is a physical primary
1512   if (!fMCevent) return kFALSE;
1513   if (fTrackLabel<0) return kFALSE;
1514   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
1515 }
1516
1517 //-----------------------------------------------------------------------
1518 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
1519 {
1520   //check if current particle is a physical primary
1521   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
1522   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
1523   if (!track) return kFALSE;
1524   TParticle* particle = track->Particle();
1525   Bool_t transported = particle->TestBit(kTransportBit);
1526   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
1527         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
1528   return (physprim && (transported || !requiretransported));
1529 }
1530
1531 //-----------------------------------------------------------------------
1532 void AliFlowTrackCuts::DefineHistograms()
1533 {
1534   //define qa histograms
1535   if (fQA) return;
1536
1537   const Int_t kNbinsP=200;
1538   Double_t binsP[kNbinsP+1];
1539   binsP[0]=0.0;
1540   for(int i=1; i<kNbinsP+1; i++)
1541   {
1542     //if(binsP[i-1]+0.05<1.01)
1543     //  binsP[i]=binsP[i-1]+0.05;
1544     //else
1545     binsP[i]=binsP[i-1]+0.05;
1546   }
1547
1548   const Int_t nBinsDCA=1000;
1549   Double_t binsDCA[nBinsDCA+1];
1550   for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
1551   //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
1552
1553   Bool_t adddirstatus = TH1::AddDirectoryStatus();
1554   TH1::AddDirectory(kFALSE);
1555   fQA=new TList(); fQA->SetOwner();
1556   fQA->SetName(Form("%s QA",GetName()));
1557   TList* before = new TList(); before->SetOwner();
1558   before->SetName("before");
1559   TList* after = new TList(); after->SetOwner();
1560   after->SetName("after");
1561   fQA->Add(before);
1562   fQA->Add(after);
1563   before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
1564   after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
1565   before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
1566   after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
1567   before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
1568   after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
1569   //primary
1570   TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
1571   TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
1572   TAxis* axis = NULL;
1573   axis = hb->GetYaxis();
1574   axis->SetBinLabel(1,"secondary");
1575   axis->SetBinLabel(2,"primary");
1576   axis = ha->GetYaxis();
1577   axis->SetBinLabel(1,"secondary");
1578   axis->SetBinLabel(2,"primary");
1579   before->Add(hb); //3
1580   after->Add(ha); //3
1581   //production process
1582   hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
1583                       -0.5, kMaxMCProcess-0.5);
1584   ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
1585                       -0.5, kMaxMCProcess-0.5);
1586   axis = hb->GetYaxis();
1587   for (Int_t i=0; i<kMaxMCProcess; i++)
1588   {
1589     axis->SetBinLabel(i+1,TMCProcessName[i]);
1590   }
1591   axis = ha->GetYaxis();
1592   for (Int_t i=0; i<kMaxMCProcess; i++)
1593   {
1594     axis->SetBinLabel(i+1,TMCProcessName[i]);
1595   }
1596   before->Add(hb); //4
1597   after->Add(ha); //4
1598   //DCA
1599   before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
1600   after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
1601   before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
1602   after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
1603   //first mother
1604   hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
1605   ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
1606   hb->GetYaxis()->SetBinLabel(1, "#gamma");
1607   ha->GetYaxis()->SetBinLabel(1, "#gamma");
1608   hb->GetYaxis()->SetBinLabel(2, "e^{+}");
1609   ha->GetYaxis()->SetBinLabel(2, "e^{+}");
1610   hb->GetYaxis()->SetBinLabel(3, "e^{-}");
1611   ha->GetYaxis()->SetBinLabel(3, "e^{-}");
1612   hb->GetYaxis()->SetBinLabel(4, "#nu");
1613   ha->GetYaxis()->SetBinLabel(4, "#nu");
1614   hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
1615   ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
1616   hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
1617   ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
1618   hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
1619   ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
1620   hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
1621   ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
1622   hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
1623   ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
1624   hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
1625   ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
1626   hb->GetYaxis()->SetBinLabel(11, "K^{+}");
1627   ha->GetYaxis()->SetBinLabel(11, "K^{+}");
1628   hb->GetYaxis()->SetBinLabel(12, "K^{-}");
1629   ha->GetYaxis()->SetBinLabel(12, "K^{-}");
1630   hb->GetYaxis()->SetBinLabel( 13, "n");
1631   ha->GetYaxis()->SetBinLabel( 13, "n");
1632   hb->GetYaxis()->SetBinLabel( 14, "p");
1633   ha->GetYaxis()->SetBinLabel( 14, "p");
1634   hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
1635   ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
1636   hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
1637   ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
1638   hb->GetYaxis()->SetBinLabel(17, "#eta");
1639   ha->GetYaxis()->SetBinLabel(17, "#eta");
1640   hb->GetYaxis()->SetBinLabel(18, "#Lambda");
1641   ha->GetYaxis()->SetBinLabel(18, "#Lambda");
1642   hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
1643   ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
1644   hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
1645   ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
1646   hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
1647   ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
1648   hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
1649   ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
1650   hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
1651   ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
1652   hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
1653   ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
1654   hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
1655   ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
1656   hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
1657   ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
1658   hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
1659   ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
1660   hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
1661   ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
1662   hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
1663   ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
1664   hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
1665   ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
1666   hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
1667   ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
1668   hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
1669   ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
1670   hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
1671   ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
1672   hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
1673   ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
1674   hb->GetYaxis()->SetBinLabel(35, "D^{+}");
1675   ha->GetYaxis()->SetBinLabel(35, "D^{+}");
1676   hb->GetYaxis()->SetBinLabel(36, "D^{-}");
1677   ha->GetYaxis()->SetBinLabel(36, "D^{-}");
1678   hb->GetYaxis()->SetBinLabel(37, "D^{0}");
1679   ha->GetYaxis()->SetBinLabel(37, "D^{0}");
1680   hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
1681   ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
1682   hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
1683   ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
1684   hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
1685   ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
1686   hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
1687   ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
1688   hb->GetYaxis()->SetBinLabel(42, "W^{+}");
1689   ha->GetYaxis()->SetBinLabel(42, "W^{+}");
1690   hb->GetYaxis()->SetBinLabel(43, "W^{-}");
1691   ha->GetYaxis()->SetBinLabel(43, "W^{-}");
1692   hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
1693   ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
1694   before->Add(hb);//7
1695   after->Add(ha);//7
1696
1697   TH1::AddDirectory(adddirstatus);
1698 }
1699
1700 //-----------------------------------------------------------------------
1701 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
1702 {
1703   //get the number of tracks in the input event according source
1704   //selection (ESD tracks, tracklets, MC particles etc.)
1705   AliESDEvent* esd=NULL;
1706   switch (fParamType)
1707   {
1708     case kSPDtracklet:
1709       esd = dynamic_cast<AliESDEvent*>(fEvent);
1710       if (!esd) return 0;
1711       return esd->GetMultiplicity()->GetNumberOfTracklets();
1712     case kMC:
1713       if (!fMCevent) return 0;
1714       return fMCevent->GetNumberOfTracks();
1715     case kPMD:
1716       esd = dynamic_cast<AliESDEvent*>(fEvent);
1717       if (!esd) return 0;
1718       return esd->GetNumberOfPmdTracks();
1719     case kV0:
1720       return fgkNumberOfV0tracks;
1721     default:
1722       if (!fEvent) return 0;
1723       return fEvent->GetNumberOfTracks();
1724   }
1725   return 0;
1726 }
1727
1728 //-----------------------------------------------------------------------
1729 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
1730 {
1731   //get the input object according the data source selection:
1732   //(esd tracks, traclets, mc particles,etc...)
1733   AliESDEvent* esd=NULL;
1734   switch (fParamType)
1735   {
1736     case kSPDtracklet:
1737       esd = dynamic_cast<AliESDEvent*>(fEvent);
1738       if (!esd) return NULL;
1739       return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
1740     case kMC:
1741       if (!fMCevent) return NULL;
1742       return fMCevent->GetTrack(i);
1743     case kPMD:
1744       esd = dynamic_cast<AliESDEvent*>(fEvent);
1745       if (!esd) return NULL;
1746       return esd->GetPmdTrack(i);
1747     case kV0:
1748       //esd = dynamic_cast<AliESDEvent*>(fEvent);
1749       //if (!esd) //contributed by G.Ortona
1750       //{
1751       //  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fEvent);
1752       //  if(!aod)return NULL;
1753       //  return aod->GetVZEROData();
1754       //}
1755       //return esd->GetVZEROData();
1756       return fEvent; // left only for compatibility
1757     default:
1758       if (!fEvent) return NULL;
1759       return fEvent->GetTrack(i);
1760   }
1761 }
1762
1763 //-----------------------------------------------------------------------
1764 void AliFlowTrackCuts::Clear(Option_t*)
1765 {
1766   //clean up
1767   fTrack=NULL;
1768   fMCevent=NULL;
1769   fMCparticle=NULL;
1770   fTrackLabel=-997;
1771   fTrackWeight=1.0;
1772   fTrackEta=0.0;
1773   fTrackPhi=0.0;
1774 }
1775
1776 //-----------------------------------------------------------------------
1777 Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
1778 {
1779   //check if passes the selected pid cut for ESDs
1780   Bool_t pass = kTRUE;
1781   switch (fPIDsource)
1782   {
1783     case kTPCpid:
1784       if (!PassesTPCpidCut(track)) pass=kFALSE;
1785       break;
1786     case kTPCdedx:
1787       if (!PassesTPCdedxCut(track)) pass=kFALSE;
1788       break;
1789     case kTOFpid:
1790       if (!PassesTOFpidCut(track)) pass=kFALSE;
1791       break;
1792     case kTOFbeta:
1793       if (!PassesTOFbetaCut(track)) pass=kFALSE;
1794       break;
1795     case kTOFbetaSimple:
1796       if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
1797       break;
1798     case kTPCbayesian:
1799       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
1800       break;
1801       // part added by F. Noferini
1802     case kTOFbayesian:
1803       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
1804       break;
1805       // end part added by F. Noferini
1806
1807       //part added by Natasha
1808     case kTPCNuclei:
1809       if (!PassesNucleiSelection(track)) pass=kFALSE;
1810       break;
1811       //end part added by Natasha
1812
1813     default:
1814       printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
1815       pass=kFALSE;
1816       break;
1817   }
1818   return pass;
1819 }
1820
1821 //-----------------------------------------------------------------------
1822 Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
1823 {
1824   //check if passes PID cut using timing in TOF
1825   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
1826                      (track->GetTOFsignal() > 12000) &&
1827                      (track->GetTOFsignal() < 100000) &&
1828                      (track->GetIntegratedLength() > 365);
1829
1830   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
1831
1832   Bool_t statusMatchingHard = TPCTOFagree(track);
1833   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
1834        return kFALSE;
1835
1836   if (!goodtrack) return kFALSE;
1837
1838   const Float_t c = 2.99792457999999984e-02;
1839   Float_t p = track->GetP();
1840   Float_t l = track->GetIntegratedLength();
1841   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
1842   Float_t timeTOF = track->GetTOFsignal()- trackT0;
1843   Float_t beta = l/timeTOF/c;
1844   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
1845   track->GetIntegratedTimes(integratedTimes);
1846   Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
1847   Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
1848   for (Int_t i=0;i<5;i++)
1849   {
1850     betaHypothesis[i] = l/integratedTimes[i]/c;
1851     s[i] = beta-betaHypothesis[i];
1852   }
1853
1854   switch (fParticleID)
1855   {
1856     case AliPID::kPion:
1857       return ( (s[2]<0.015) && (s[2]>-0.015) &&
1858                (s[3]>0.025) &&
1859                (s[4]>0.03) );
1860     case AliPID::kKaon:
1861       return ( (s[3]<0.015) && (s[3]>-0.015) &&
1862                (s[2]<-0.03) &&
1863                (s[4]>0.03) );
1864     case AliPID::kProton:
1865       return ( (s[4]<0.015) && (s[4]>-0.015) &&
1866                (s[3]<-0.025) &&
1867                (s[2]<-0.025) );
1868     default:
1869       return kFALSE;
1870   }
1871   return kFALSE;
1872 }
1873
1874 //-----------------------------------------------------------------------
1875 Float_t AliFlowTrackCuts::GetBeta(const AliESDtrack* track)
1876 {
1877   //get beta
1878   const Float_t c = 2.99792457999999984e-02;
1879   Float_t p = track->GetP();
1880   Float_t l = track->GetIntegratedLength();
1881   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
1882   Float_t timeTOF = track->GetTOFsignal()- trackT0;
1883   return l/timeTOF/c;
1884 }
1885
1886 //-----------------------------------------------------------------------
1887 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
1888 {
1889   //check if track passes pid selection with an asymmetric TOF beta cut
1890   if (!fTOFpidCuts)
1891   {
1892     //printf("no TOFpidCuts\n");
1893     return kFALSE;
1894   }
1895
1896   //check if passes PID cut using timing in TOF
1897   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
1898                      (track->GetTOFsignal() > 12000) &&
1899                      (track->GetTOFsignal() < 100000) &&
1900                      (track->GetIntegratedLength() > 365);
1901
1902   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
1903
1904   Bool_t statusMatchingHard = TPCTOFagree(track);
1905   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
1906        return kFALSE;
1907
1908   if (!goodtrack) return kFALSE;
1909
1910   Float_t beta = GetBeta(track);
1911
1912   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
1913   track->GetIntegratedTimes(integratedTimes);
1914
1915   //construct the pid index because it's not AliPID::EParticleType
1916   Int_t pid = 0;
1917   switch (fParticleID)
1918   {
1919     case AliPID::kPion:
1920       pid=2;
1921       break;
1922     case AliPID::kKaon:
1923       pid=3;
1924       break;
1925     case AliPID::kProton:
1926       pid=4;
1927       break;
1928     default:
1929       return kFALSE;
1930   }
1931
1932   //signal to cut on
1933   const Float_t c = 2.99792457999999984e-02;
1934   Float_t l = track->GetIntegratedLength();
1935   Float_t p = track->GetP();
1936   Float_t betahypothesis = l/integratedTimes[pid]/c;
1937   Float_t betadiff = beta-betahypothesis;
1938
1939   Float_t* arr = fTOFpidCuts->GetMatrixArray();
1940   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
1941   if (col<0) return kFALSE;
1942   Float_t min = (*fTOFpidCuts)(1,col);
1943   Float_t max = (*fTOFpidCuts)(2,col);
1944
1945   Bool_t pass = (betadiff>min && betadiff<max);
1946
1947   return pass;
1948 }
1949
1950 //-----------------------------------------------------------------------
1951 Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
1952 {
1953   //check if passes PID cut using default TOF pid
1954   Double_t pidTOF[AliPID::kSPECIES];
1955   track->GetTOFpid(pidTOF);
1956   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
1957   return kFALSE;
1958 }
1959
1960 //-----------------------------------------------------------------------
1961 Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
1962 {
1963   //check if passes PID cut using default TPC pid
1964   Double_t pidTPC[AliPID::kSPECIES];
1965   track->GetTPCpid(pidTPC);
1966   Double_t probablity = 0.;
1967   switch (fParticleID)
1968   {
1969     case AliPID::kPion:
1970       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
1971       break;
1972     default:
1973       probablity = pidTPC[fParticleID];
1974   }
1975   if (probablity >= fParticleProbability) return kTRUE;
1976   return kFALSE;
1977 }
1978
1979 //-----------------------------------------------------------------------
1980 Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
1981 {
1982   //get TPC dedx
1983   return track->GetTPCsignal();
1984 }
1985
1986 //-----------------------------------------------------------------------
1987 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
1988 {
1989   //check if passes PID cut using dedx signal in the TPC
1990   if (!fTPCpidCuts)
1991   {
1992     //printf("no TPCpidCuts\n");
1993     return kFALSE;
1994   }
1995
1996   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
1997   if (!tpcparam) return kFALSE;
1998   Double_t p = tpcparam->GetP();
1999   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
2000   Float_t sigTPC = track->GetTPCsignal();
2001   Float_t s = (sigTPC-sigExp)/sigExp;
2002
2003   Float_t* arr = fTPCpidCuts->GetMatrixArray();
2004   Int_t arrSize = fTPCpidCuts->GetNcols();
2005   Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
2006   if (col<0) return kFALSE;
2007   Float_t min = (*fTPCpidCuts)(1,col);
2008   Float_t max = (*fTPCpidCuts)(2,col);
2009
2010   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
2011   return (s>min && s<max);
2012 }
2013
2014 //-----------------------------------------------------------------------
2015 void AliFlowTrackCuts::InitPIDcuts()
2016 {
2017   //init matrices with PID cuts
2018   TMatrixF* t = NULL;
2019   if (!fTPCpidCuts)
2020   {
2021     if (fParticleID==AliPID::kPion)
2022     {
2023       t = new TMatrixF(3,15);
2024       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.0;
2025       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.1;
2026       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.2;
2027       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.2;
2028       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
2029       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
2030       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
2031       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
2032       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
2033       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.4;  (*t)(2,9)  =  0.05;
2034       (*t)(0,10)  = 0.70;  (*t)(1,10)  = -0.4;  (*t)(2,10)  =     0;
2035       (*t)(0,11)  = 0.75;  (*t)(1,11)  = -0.4;  (*t)(2,11)  =     0;
2036       (*t)(0,12)  = 0.80;  (*t)(1,12)  = -0.4;  (*t)(2,12)  = -0.05;
2037       (*t)(0,13)  = 0.85;  (*t)(1,13)  = -0.4;  (*t)(2,13)  = -0.1;
2038       (*t)(0,14)  = 0.90;  (*t)(1,14)  = 0;     (*t)(2,14)  =     0;
2039     }
2040     else
2041     if (fParticleID==AliPID::kKaon)
2042     {
2043       t = new TMatrixF(3,12);
2044       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.2;
2045       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  = 0.2;
2046       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.2;  (*t)(2,2)  = 0.2;
2047       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.2;  (*t)(2,3)  = 0.2;
2048       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.2;
2049       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.2;
2050       (*t)(0,6)  = 0.50;  (*t)(1,6)  =-0.05;  (*t)(2,6)  = 0.2;
2051       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.1;  (*t)(2,7)  = 0.1;
2052       (*t)(0,8)  = 0.60;  (*t)(1,8)  =-0.05;  (*t)(2,8)  = 0.1;
2053       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  = 0.15;
2054       (*t)(0,10)  = 0.70;  (*t)(1,10)  = 0.05;  (*t)(2,10)  = 0.2;
2055       (*t)(0,11)  = 0.75;  (*t)(1,11)  =    0;  (*t)(2,11)  = 0;
2056     }
2057     else
2058     if (fParticleID==AliPID::kProton)
2059     {
2060       t = new TMatrixF(3,9);
2061       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.1;  (*t)(2,0)  =  0.1;
2062       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  =  0.2;
2063       (*t)(0,2)  = 0.80;  (*t)(1,2)  = -0.1;  (*t)(2,2)  =  0.2;
2064       (*t)(0,3)  = 0.85;  (*t)(1,3)  =-0.05;  (*t)(2,3)  =  0.2;
2065       (*t)(0,4)  = 0.90;  (*t)(1,4)  =-0.05;  (*t)(2,4)  = 0.25;
2066       (*t)(0,5)  = 0.95;  (*t)(1,5)  =-0.05;  (*t)(2,5)  = 0.25;
2067       (*t)(0,6)  = 1.00;  (*t)(1,6)  = -0.1;  (*t)(2,6)  = 0.25;
2068       (*t)(0,7)  = 1.10;  (*t)(1,7)  =-0.05;  (*t)(2,7)  =  0.3;
2069       (*t)(0,8) = 1.20;   (*t)(1,8)  =    0;  (*t)(2,8) =    0;
2070     }
2071     delete fTPCpidCuts;
2072     fTPCpidCuts=t;
2073   }
2074   t = NULL;
2075   if (!fTOFpidCuts)
2076   {
2077     if (fParticleID==AliPID::kPion)
2078     {
2079       //TOF pions, 0.9 purity
2080       t = new TMatrixF(3,61);
2081       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2082       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2083       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2084       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2085       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.030;  (*t)(2,4)  =   0.030;
2086       (*t)(0,5)  = 0.250;  (*t)(1,5)  = -0.036;  (*t)(2,5)  =   0.032;
2087       (*t)(0,6)  = 0.300;  (*t)(1,6)  = -0.038;  (*t)(2,6)  =   0.032;
2088       (*t)(0,7)  = 0.350;  (*t)(1,7)  = -0.034;  (*t)(2,7)  =   0.032;
2089       (*t)(0,8)  = 0.400;  (*t)(1,8)  = -0.032;  (*t)(2,8)  =   0.020;
2090       (*t)(0,9)  = 0.450;  (*t)(1,9)  = -0.030;  (*t)(2,9)  =   0.020;
2091       (*t)(0,10)  = 0.500;  (*t)(1,10)  = -0.030;  (*t)(2,10)  =   0.020;
2092       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.030;  (*t)(2,11)  =   0.020;
2093       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.030;  (*t)(2,12)  =   0.020;
2094       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.030;  (*t)(2,13)  =   0.020;
2095       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.030;  (*t)(2,14)  =   0.020;
2096       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.030;  (*t)(2,15)  =   0.020;
2097       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.030;  (*t)(2,16)  =   0.020;
2098       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.030;  (*t)(2,17)  =   0.020;
2099       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.030;  (*t)(2,18)  =   0.020;
2100       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.028;  (*t)(2,19)  =   0.028;
2101       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.028;  (*t)(2,20)  =   0.028;
2102       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.028;  (*t)(2,21)  =   0.028;
2103       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.028;
2104       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.028;
2105       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.020;  (*t)(2,24)  =   0.028;
2106       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.018;  (*t)(2,25)  =   0.028;
2107       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.016;  (*t)(2,26)  =   0.028;
2108       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.014;  (*t)(2,27)  =   0.028;
2109       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.012;  (*t)(2,28)  =   0.026;
2110       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.010;  (*t)(2,29)  =   0.026;
2111       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.008;  (*t)(2,30)  =   0.026;
2112       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.008;  (*t)(2,31)  =   0.024;
2113       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.006;  (*t)(2,32)  =   0.024;
2114       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.004;  (*t)(2,33)  =   0.024;
2115       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.004;  (*t)(2,34)  =   0.024;
2116       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.002;  (*t)(2,35)  =   0.024;
2117       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.002;  (*t)(2,36)  =   0.024;
2118       (*t)(0,37)  = 2.700;  (*t)(1,37)  = 0.000;  (*t)(2,37)  =   0.024;
2119       (*t)(0,38)  = 2.800;  (*t)(1,38)  = 0.000;  (*t)(2,38)  =   0.026;
2120       (*t)(0,39)  = 2.900;  (*t)(1,39)  = 0.000;  (*t)(2,39)  =   0.024;
2121       (*t)(0,40)  = 3.000;  (*t)(1,40)  = 0.002;  (*t)(2,40)  =   0.026;
2122       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.002;  (*t)(2,41)  =   0.026;
2123       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.002;  (*t)(2,42)  =   0.026;
2124       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.002;  (*t)(2,43)  =   0.026;
2125       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.002;  (*t)(2,44)  =   0.026;
2126       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.002;  (*t)(2,45)  =   0.026;
2127       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.002;  (*t)(2,46)  =   0.026;
2128       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.002;  (*t)(2,47)  =   0.026;
2129       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.002;  (*t)(2,48)  =   0.026;
2130       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.004;  (*t)(2,49)  =   0.024;
2131       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.004;  (*t)(2,50)  =   0.026;
2132       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.004;  (*t)(2,51)  =   0.026;
2133       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.004;  (*t)(2,52)  =   0.024;
2134       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.006;  (*t)(2,53)  =   0.024;
2135       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2136       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2137       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2138       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2139       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2140       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2141       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2142     }
2143     else
2144     if (fParticleID==AliPID::kProton)
2145     {
2146       //TOF protons, 0.9 purity
2147       t = new TMatrixF(3,61);
2148       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2149       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2150       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2151       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2152       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.07;  (*t)(2,4)  =   0.07;
2153       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.07;  (*t)(2,5)  =   0.07;
2154       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.07;  (*t)(2,6)  =   0.07;
2155       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.07;  (*t)(2,7)  =   0.07;
2156       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.07;  (*t)(2,8)  =   0.07;
2157       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.07;  (*t)(2,9)  =   0.07;
2158       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.07;  (*t)(2,10)  =   0.07;
2159       (*t)(0,11)  = 0.200;  (*t)(1,11)  = -0.07;  (*t)(2,11)  =   0.07;
2160       (*t)(0,12)  = 0.200;  (*t)(1,12)  = -0.07;  (*t)(2,12)  =   0.07;
2161       (*t)(0,13)  = 0.200;  (*t)(1,13)  = -0.07;  (*t)(2,13)  =   0.07;
2162       (*t)(0,14)  = 0.200;  (*t)(1,14)  = -0.07;  (*t)(2,14)  =   0.07;
2163       (*t)(0,15)  = 0.200;  (*t)(1,15)  = -0.07;  (*t)(2,15)  =   0.07;
2164       (*t)(0,16)  = 0.200;  (*t)(1,16)  = -0.07;  (*t)(2,16)  =   0.07;
2165       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.070;  (*t)(2,17)  =   0.070;
2166       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.072;  (*t)(2,18)  =   0.072;
2167       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.072;  (*t)(2,19)  =   0.072;
2168       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.074;  (*t)(2,20)  =   0.074;
2169       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.032;  (*t)(2,21)  =   0.032;
2170       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.026;
2171       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.026;  (*t)(2,23)  =   0.026;
2172       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.024;  (*t)(2,24)  =   0.024;
2173       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.024;  (*t)(2,25)  =   0.024;
2174       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.026;  (*t)(2,26)  =   0.026;
2175       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.026;  (*t)(2,27)  =   0.026;
2176       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.026;  (*t)(2,28)  =   0.026;
2177       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.026;  (*t)(2,29)  =   0.026;
2178       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.026;  (*t)(2,30)  =   0.026;
2179       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.026;
2180       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.026;  (*t)(2,32)  =   0.024;
2181       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.028;  (*t)(2,33)  =   0.022;
2182       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.028;  (*t)(2,34)  =   0.020;
2183       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.028;  (*t)(2,35)  =   0.018;
2184       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.028;  (*t)(2,36)  =   0.016;
2185       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.028;  (*t)(2,37)  =   0.016;
2186       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.030;  (*t)(2,38)  =   0.014;
2187       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.030;  (*t)(2,39)  =   0.012;
2188       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.030;  (*t)(2,40)  =   0.012;
2189       (*t)(0,41)  = 3.100;  (*t)(1,41)  = -0.030;  (*t)(2,41)  =   0.010;
2190       (*t)(0,42)  = 3.200;  (*t)(1,42)  = -0.030;  (*t)(2,42)  =   0.010;
2191       (*t)(0,43)  = 3.300;  (*t)(1,43)  = -0.030;  (*t)(2,43)  =   0.010;
2192       (*t)(0,44)  = 3.400;  (*t)(1,44)  = -0.030;  (*t)(2,44)  =   0.008;
2193       (*t)(0,45)  = 3.500;  (*t)(1,45)  = -0.030;  (*t)(2,45)  =   0.008;
2194       (*t)(0,46)  = 3.600;  (*t)(1,46)  = -0.030;  (*t)(2,46)  =   0.008;
2195       (*t)(0,47)  = 3.700;  (*t)(1,47)  = -0.030;  (*t)(2,47)  =   0.006;
2196       (*t)(0,48)  = 3.800;  (*t)(1,48)  = -0.030;  (*t)(2,48)  =   0.006;
2197       (*t)(0,49)  = 3.900;  (*t)(1,49)  = -0.030;  (*t)(2,49)  =   0.006;
2198       (*t)(0,50)  = 4.000;  (*t)(1,50)  = -0.028;  (*t)(2,50)  =   0.004;
2199       (*t)(0,51)  = 4.100;  (*t)(1,51)  = -0.030;  (*t)(2,51)  =   0.004;
2200       (*t)(0,52)  = 4.200;  (*t)(1,52)  = -0.030;  (*t)(2,52)  =   0.004;
2201       (*t)(0,53)  = 4.300;  (*t)(1,53)  = -0.028;  (*t)(2,53)  =   0.002;
2202       (*t)(0,54)  = 4.400;  (*t)(1,54)  = -0.030;  (*t)(2,54)  =   0.002;
2203       (*t)(0,55)  = 4.500;  (*t)(1,55)  = -0.028;  (*t)(2,55)  =   0.002;
2204       (*t)(0,56)  = 4.600;  (*t)(1,56)  = -0.028;  (*t)(2,56)  =   0.002;
2205       (*t)(0,57)  = 4.700;  (*t)(1,57)  = -0.028;  (*t)(2,57)  =   0.000;
2206       (*t)(0,58)  = 4.800;  (*t)(1,58)  = -0.028;  (*t)(2,58)  =   0.002;
2207       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2208       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2209     }
2210     else
2211     if (fParticleID==AliPID::kKaon)
2212     {
2213       //TOF kaons, 0.9 purity
2214       t = new TMatrixF(3,61);
2215       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2216       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2217       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2218       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2219       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.05;  (*t)(2,4)  =   0.05;
2220       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.05;  (*t)(2,5)  =   0.05;
2221       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.05;  (*t)(2,6)  =   0.05;
2222       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =   0.05;
2223       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.05;  (*t)(2,8)  =   0.05;
2224       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.05;  (*t)(2,9)  =   0.05;
2225       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.05;  (*t)(2,10)  =   0.05;
2226       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.026;  (*t)(2,11)  =   0.026;
2227       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.026;  (*t)(2,12)  =   0.026;
2228       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.026;  (*t)(2,13)  =   0.026;
2229       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.026;  (*t)(2,14)  =   0.026;
2230       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.026;  (*t)(2,15)  =   0.026;
2231       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.026;  (*t)(2,16)  =   0.026;
2232       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.024;  (*t)(2,17)  =   0.024;
2233       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.024;  (*t)(2,18)  =   0.024;
2234       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.024;  (*t)(2,19)  =   0.024;
2235       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.024;  (*t)(2,20)  =   0.024;
2236       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.024;  (*t)(2,21)  =   0.024;
2237       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.024;  (*t)(2,22)  =   0.022;
2238       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.020;
2239       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.026;  (*t)(2,24)  =   0.016;
2240       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.028;  (*t)(2,25)  =   0.014;
2241       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.028;  (*t)(2,26)  =   0.012;
2242       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.028;  (*t)(2,27)  =   0.010;
2243       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.028;  (*t)(2,28)  =   0.010;
2244       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.028;  (*t)(2,29)  =   0.008;
2245       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.028;  (*t)(2,30)  =   0.006;
2246       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.006;
2247       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.024;  (*t)(2,32)  =   0.004;
2248       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.020;  (*t)(2,33)  =   0.002;
2249       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.020;  (*t)(2,34)  =   0.002;
2250       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.018;  (*t)(2,35)  =   0.000;
2251       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.016;  (*t)(2,36)  =   0.000;
2252       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.014;  (*t)(2,37)  =   -0.002;
2253       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.014;  (*t)(2,38)  =   -0.004;
2254       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.012;  (*t)(2,39)  =   -0.004;
2255       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.010;  (*t)(2,40)  =   -0.006;
2256       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.000;  (*t)(2,41)  =   0.000;
2257       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.000;  (*t)(2,42)  =   0.000;
2258       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.000;  (*t)(2,43)  =   0.000;
2259       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.000;  (*t)(2,44)  =   0.000;
2260       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.000;  (*t)(2,45)  =   0.000;
2261       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.000;  (*t)(2,46)  =   0.000;
2262       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.000;  (*t)(2,47)  =   0.000;
2263       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.000;  (*t)(2,48)  =   0.000;
2264       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.000;  (*t)(2,49)  =   0.000;
2265       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.000;  (*t)(2,50)  =   0.000;
2266       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.000;  (*t)(2,51)  =   0.000;
2267       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.000;  (*t)(2,52)  =   0.000;
2268       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.000;  (*t)(2,53)  =   0.000;
2269       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2270       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2271       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2272       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2273       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2274       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2275       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2276     }
2277     delete fTOFpidCuts;
2278     fTOFpidCuts=t;
2279   }
2280 }
2281
2282 //-----------------------------------------------------------------------
2283 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
2284 {
2285   //cut on TPC bayesian pid
2286   //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
2287
2288   //Bool_t statusMatchingHard = TPCTOFagree(track);
2289   //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2290   //     return kFALSE;
2291   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
2292   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
2293   Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
2294
2295   if(! kTPC) return kFALSE;
2296
2297   Bool_t statusMatchingHard = 1;
2298   Float_t mismProb = 0;
2299   if(kTOF){
2300     statusMatchingHard = TPCTOFagree(track);
2301     mismProb = fBayesianResponse->GetTOFMismProb();
2302   }
2303   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2304        return kFALSE;
2305
2306   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
2307
2308   fProbBayes = 0.0;
2309
2310   switch (fParticleID)
2311   {
2312     case AliPID::kPion:
2313       fProbBayes = probabilities[2];
2314       break;
2315     case AliPID::kKaon:
2316       fProbBayes = probabilities[3];
2317      break;
2318     case AliPID::kProton:
2319       fProbBayes = probabilities[4];
2320       break;
2321     case AliPID::kElectron:
2322        fProbBayes = probabilities[0];
2323      break;
2324     case AliPID::kMuon:
2325        fProbBayes = probabilities[1];
2326      break;
2327     case AliPID::kDeuteron:
2328        fProbBayes = probabilities[5];
2329      break;
2330     case AliPID::kTriton:
2331        fProbBayes = probabilities[6];
2332      break;
2333     case AliPID::kHe3:
2334        fProbBayes = probabilities[7];
2335      break;
2336     default:
2337       return kFALSE;
2338   }
2339
2340   if(fProbBayes > fParticleProbability && mismProb < 0.5)
2341     {
2342       if(!fCutCharge)
2343         return kTRUE;
2344       else if (fCutCharge && fCharge * track->GetSign() > 0)
2345         return kTRUE;
2346     }
2347   return kFALSE;
2348 }
2349
2350 //-----------------------------------------------------------------------
2351 // part added by F. Noferini (some methods)
2352 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
2353 {
2354   //check is track passes bayesian combined TOF+TPC pid cut
2355   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
2356                      (track->GetStatus() & AliESDtrack::kTIME) &&
2357                      (track->GetTOFsignal() > 12000) &&
2358                      (track->GetTOFsignal() < 100000) &&
2359                      (track->GetIntegratedLength() > 365);
2360
2361   if (! goodtrack)
2362        return kFALSE;
2363
2364   if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
2365
2366   Bool_t statusMatchingHard = TPCTOFagree(track);
2367   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2368        return kFALSE;
2369
2370   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
2371   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
2372
2373   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2374
2375   fProbBayes = 0.0;
2376
2377   switch (fParticleID)
2378   {
2379     case AliPID::kPion:
2380       fProbBayes = probabilities[2];
2381       break;
2382     case AliPID::kKaon:
2383        fProbBayes = probabilities[3];
2384      break;
2385     case AliPID::kProton:
2386        fProbBayes = probabilities[4];
2387       break;
2388     case AliPID::kElectron:
2389        fProbBayes = probabilities[0];
2390      break;
2391     case AliPID::kMuon:
2392        fProbBayes = probabilities[1];
2393      break;
2394     case AliPID::kDeuteron:
2395        fProbBayes = probabilities[5];
2396      break;
2397     case AliPID::kTriton:
2398        fProbBayes = probabilities[6];
2399      break;
2400     case AliPID::kHe3:
2401        fProbBayes = probabilities[7];
2402      break;
2403    default:
2404       return kFALSE;
2405   }
2406
2407   //  printf("pt = %f -- all prob = [%4.2f,%4.2f,%4.2f,%4.2f,%4.2f] -- prob = %f\n",track->Pt(),fProbBayes[0],fProbBayes[1],fProbBayes[2],fProbBayes[3],fProbBayes[4],prob);
2408   if(fProbBayes > fParticleProbability && mismProb < 0.5){
2409     if(!fCutCharge)
2410       return kTRUE;
2411     else if (fCutCharge && fCharge * track->GetSign() > 0)
2412       return kTRUE;
2413   }
2414   return kFALSE;
2415 }
2416
2417
2418 //-----------------------------------------------------------------------
2419  // part added by Natasha
2420 Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
2421 {
2422   //pid selection for heavy nuclei
2423   Bool_t select=kFALSE;
2424
2425   //if (!track) continue;
2426
2427   if (!track->GetInnerParam())
2428     return kFALSE;    //break;
2429
2430   const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
2431
2432   Double_t ptotTPC = tpcTrack->GetP();
2433   Double_t sigTPC = track->GetTPCsignal();
2434   Double_t dEdxBBA = 0.;
2435   Double_t dSigma = 0.;
2436
2437   switch (fParticleID)
2438   {
2439     case AliPID::kDeuteron:
2440       //pid=10;
2441       dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
2442           4.60e+00,
2443           8.9684e+00,
2444           1.640e-05,
2445           2.35e+00,
2446           2.35e+00);
2447       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
2448
2449       if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) )  )
2450       {select=kTRUE;}
2451       break;
2452
2453     case AliPID::kTriton:
2454       //pid=11;
2455       select=kFALSE;
2456       break;
2457
2458     case AliPID::kHe3:
2459       //pid=12;
2460       // ----- Pass 2 -------
2461       dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
2462           1.74962,
2463           27.4992,
2464           4.00313e-15,
2465           2.42485,
2466           8.31768);
2467       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
2468       if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
2469       {select=kTRUE;}
2470       break;
2471
2472     case AliPID::kAlpha:
2473       //pid=13;
2474       select=kFALSE;
2475       break;
2476
2477     default:
2478       return kFALSE;
2479   }
2480
2481   return select;
2482 }
2483 // end part added by Natasha
2484 //-----------------------------------------------------------------------
2485 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
2486  //set priors for the bayesian pid selection
2487   fCurrCentr = centrCur;
2488
2489   fBinLimitPID[0] = 0.300000;
2490   fBinLimitPID[1] = 0.400000;
2491   fBinLimitPID[2] = 0.500000;
2492   fBinLimitPID[3] = 0.600000;
2493   fBinLimitPID[4] = 0.700000;
2494   fBinLimitPID[5] = 0.800000;
2495   fBinLimitPID[6] = 0.900000;
2496   fBinLimitPID[7] = 1.000000;
2497   fBinLimitPID[8] = 1.200000;
2498   fBinLimitPID[9] = 1.400000;
2499   fBinLimitPID[10] = 1.600000;
2500   fBinLimitPID[11] = 1.800000;
2501   fBinLimitPID[12] = 2.000000;
2502   fBinLimitPID[13] = 2.200000;
2503   fBinLimitPID[14] = 2.400000;
2504   fBinLimitPID[15] = 2.600000;
2505   fBinLimitPID[16] = 2.800000;
2506   fBinLimitPID[17] = 3.000000;
2507
2508   // 0-10%
2509   if(centrCur < 10){
2510       fC[0][0] = 0.005;
2511       fC[0][1] = 0.005;
2512       fC[0][2] = 1.0000;
2513       fC[0][3] = 0.0010;
2514       fC[0][4] = 0.0010;
2515
2516       fC[1][0] = 0.005;
2517       fC[1][1] = 0.005;
2518       fC[1][2] = 1.0000;
2519       fC[1][3] = 0.0168;
2520       fC[1][4] = 0.0122;
2521
2522       fC[2][0] = 0.005;
2523       fC[2][1] = 0.005;
2524       fC[2][2] = 1.0000;
2525       fC[2][3] = 0.0272;
2526       fC[2][4] = 0.0070;
2527
2528       fC[3][0] = 0.005;
2529       fC[3][1] = 0.005;
2530       fC[3][2] = 1.0000;
2531       fC[3][3] = 0.0562;
2532       fC[3][4] = 0.0258;
2533
2534       fC[4][0] = 0.005;
2535       fC[4][1] = 0.005;
2536       fC[4][2] = 1.0000;
2537       fC[4][3] = 0.0861;
2538       fC[4][4] = 0.0496;
2539
2540       fC[5][0] = 0.005;
2541       fC[5][1] = 0.005;
2542       fC[5][2] = 1.0000;
2543       fC[5][3] = 0.1168;
2544       fC[5][4] = 0.0740;
2545
2546       fC[6][0] = 0.005;
2547       fC[6][1] = 0.005;
2548       fC[6][2] = 1.0000;
2549       fC[6][3] = 0.1476;
2550       fC[6][4] = 0.0998;
2551
2552       fC[7][0] = 0.005;
2553       fC[7][1] = 0.005;
2554       fC[7][2] = 1.0000;
2555       fC[7][3] = 0.1810;
2556       fC[7][4] = 0.1296;
2557
2558       fC[8][0] = 0.005;
2559       fC[8][1] = 0.005;
2560       fC[8][2] = 1.0000;
2561       fC[8][3] = 0.2240;
2562       fC[8][4] = 0.1827;
2563
2564       fC[9][0] = 0.005;
2565       fC[9][1] = 0.005;
2566       fC[9][2] = 1.0000;
2567       fC[9][3] = 0.2812;
2568       fC[9][4] = 0.2699;
2569
2570       fC[10][0] = 0.005;
2571       fC[10][1] = 0.005;
2572       fC[10][2] = 1.0000;
2573       fC[10][3] = 0.3328;
2574       fC[10][4] = 0.3714;
2575
2576       fC[11][0] = 0.005;
2577       fC[11][1] = 0.005;
2578       fC[11][2] = 1.0000;
2579       fC[11][3] = 0.3780;
2580       fC[11][4] = 0.4810;
2581
2582       fC[12][0] = 0.005;
2583       fC[12][1] = 0.005;
2584       fC[12][2] = 1.0000;
2585       fC[12][3] = 0.4125;
2586       fC[12][4] = 0.5771;
2587
2588       fC[13][0] = 0.005;
2589       fC[13][1] = 0.005;
2590       fC[13][2] = 1.0000;
2591       fC[13][3] = 0.4486;
2592       fC[13][4] = 0.6799;
2593
2594       fC[14][0] = 0.005;
2595       fC[14][1] = 0.005;
2596       fC[14][2] = 1.0000;
2597       fC[14][3] = 0.4840;
2598       fC[14][4] = 0.7668;
2599
2600       fC[15][0] = 0.005;
2601       fC[15][1] = 0.005;
2602       fC[15][2] = 1.0000;
2603       fC[15][3] = 0.4971;
2604       fC[15][4] = 0.8288;
2605
2606       fC[16][0] = 0.005;
2607       fC[16][1] = 0.005;
2608       fC[16][2] = 1.0000;
2609       fC[16][3] = 0.4956;
2610       fC[16][4] = 0.8653;
2611
2612       fC[17][0] = 0.005;
2613       fC[17][1] = 0.005;
2614       fC[17][2] = 1.0000;
2615       fC[17][3] = 0.5173;
2616       fC[17][4] = 0.9059;
2617   }
2618   // 10-20%
2619   else if(centrCur < 20){
2620      fC[0][0] = 0.005;
2621       fC[0][1] = 0.005;
2622       fC[0][2] = 1.0000;
2623       fC[0][3] = 0.0010;
2624       fC[0][4] = 0.0010;
2625
2626       fC[1][0] = 0.005;
2627       fC[1][1] = 0.005;
2628       fC[1][2] = 1.0000;
2629       fC[1][3] = 0.0132;
2630       fC[1][4] = 0.0088;
2631
2632       fC[2][0] = 0.005;
2633       fC[2][1] = 0.005;
2634       fC[2][2] = 1.0000;
2635       fC[2][3] = 0.0283;
2636       fC[2][4] = 0.0068;
2637
2638       fC[3][0] = 0.005;
2639       fC[3][1] = 0.005;
2640       fC[3][2] = 1.0000;
2641       fC[3][3] = 0.0577;
2642       fC[3][4] = 0.0279;
2643
2644       fC[4][0] = 0.005;
2645       fC[4][1] = 0.005;
2646       fC[4][2] = 1.0000;
2647       fC[4][3] = 0.0884;
2648       fC[4][4] = 0.0534;
2649
2650       fC[5][0] = 0.005;
2651       fC[5][1] = 0.005;
2652       fC[5][2] = 1.0000;
2653       fC[5][3] = 0.1179;
2654       fC[5][4] = 0.0794;
2655
2656       fC[6][0] = 0.005;
2657       fC[6][1] = 0.005;
2658       fC[6][2] = 1.0000;
2659       fC[6][3] = 0.1480;
2660       fC[6][4] = 0.1058;
2661
2662       fC[7][0] = 0.005;
2663       fC[7][1] = 0.005;
2664       fC[7][2] = 1.0000;
2665       fC[7][3] = 0.1807;
2666       fC[7][4] = 0.1366;
2667
2668       fC[8][0] = 0.005;
2669       fC[8][1] = 0.005;
2670       fC[8][2] = 1.0000;
2671       fC[8][3] = 0.2219;
2672       fC[8][4] = 0.1891;
2673
2674       fC[9][0] = 0.005;
2675       fC[9][1] = 0.005;
2676       fC[9][2] = 1.0000;
2677       fC[9][3] = 0.2804;
2678       fC[9][4] = 0.2730;
2679
2680       fC[10][0] = 0.005;
2681       fC[10][1] = 0.005;
2682       fC[10][2] = 1.0000;
2683       fC[10][3] = 0.3283;
2684       fC[10][4] = 0.3660;
2685
2686       fC[11][0] = 0.005;
2687       fC[11][1] = 0.005;
2688       fC[11][2] = 1.0000;
2689       fC[11][3] = 0.3710;
2690       fC[11][4] = 0.4647;
2691
2692       fC[12][0] = 0.005;
2693       fC[12][1] = 0.005;
2694       fC[12][2] = 1.0000;
2695       fC[12][3] = 0.4093;
2696       fC[12][4] = 0.5566;
2697
2698       fC[13][0] = 0.005;
2699       fC[13][1] = 0.005;
2700       fC[13][2] = 1.0000;
2701       fC[13][3] = 0.4302;
2702       fC[13][4] = 0.6410;
2703
2704       fC[14][0] = 0.005;
2705       fC[14][1] = 0.005;
2706       fC[14][2] = 1.0000;
2707       fC[14][3] = 0.4649;
2708       fC[14][4] = 0.7055;
2709
2710       fC[15][0] = 0.005;
2711       fC[15][1] = 0.005;
2712       fC[15][2] = 1.0000;
2713       fC[15][3] = 0.4523;
2714       fC[15][4] = 0.7440;
2715
2716       fC[16][0] = 0.005;
2717       fC[16][1] = 0.005;
2718       fC[16][2] = 1.0000;
2719       fC[16][3] = 0.4591;
2720       fC[16][4] = 0.7799;
2721
2722       fC[17][0] = 0.005;
2723       fC[17][1] = 0.005;
2724       fC[17][2] = 1.0000;
2725       fC[17][3] = 0.4804;
2726       fC[17][4] = 0.8218;
2727   }
2728   // 20-30%
2729   else if(centrCur < 30){
2730      fC[0][0] = 0.005;
2731       fC[0][1] = 0.005;
2732       fC[0][2] = 1.0000;
2733       fC[0][3] = 0.0010;
2734       fC[0][4] = 0.0010;
2735
2736       fC[1][0] = 0.005;
2737       fC[1][1] = 0.005;
2738       fC[1][2] = 1.0000;
2739       fC[1][3] = 0.0102;
2740       fC[1][4] = 0.0064;
2741
2742       fC[2][0] = 0.005;
2743       fC[2][1] = 0.005;
2744       fC[2][2] = 1.0000;
2745       fC[2][3] = 0.0292;
2746       fC[2][4] = 0.0066;
2747
2748       fC[3][0] = 0.005;
2749       fC[3][1] = 0.005;
2750       fC[3][2] = 1.0000;
2751       fC[3][3] = 0.0597;
2752       fC[3][4] = 0.0296;
2753
2754       fC[4][0] = 0.005;
2755       fC[4][1] = 0.005;
2756       fC[4][2] = 1.0000;
2757       fC[4][3] = 0.0900;
2758       fC[4][4] = 0.0589;
2759
2760       fC[5][0] = 0.005;
2761       fC[5][1] = 0.005;
2762       fC[5][2] = 1.0000;
2763       fC[5][3] = 0.1199;
2764       fC[5][4] = 0.0859;
2765
2766       fC[6][0] = 0.005;
2767       fC[6][1] = 0.005;
2768       fC[6][2] = 1.0000;
2769       fC[6][3] = 0.1505;
2770       fC[6][4] = 0.1141;
2771
2772       fC[7][0] = 0.005;
2773       fC[7][1] = 0.005;
2774       fC[7][2] = 1.0000;
2775       fC[7][3] = 0.1805;
2776       fC[7][4] = 0.1454;
2777
2778       fC[8][0] = 0.005;
2779       fC[8][1] = 0.005;
2780       fC[8][2] = 1.0000;
2781       fC[8][3] = 0.2221;
2782       fC[8][4] = 0.2004;
2783
2784       fC[9][0] = 0.005;
2785       fC[9][1] = 0.005;
2786       fC[9][2] = 1.0000;
2787       fC[9][3] = 0.2796;
2788       fC[9][4] = 0.2838;
2789
2790       fC[10][0] = 0.005;
2791       fC[10][1] = 0.005;
2792       fC[10][2] = 1.0000;
2793       fC[10][3] = 0.3271;
2794       fC[10][4] = 0.3682;
2795
2796       fC[11][0] = 0.005;
2797       fC[11][1] = 0.005;
2798       fC[11][2] = 1.0000;
2799       fC[11][3] = 0.3648;
2800       fC[11][4] = 0.4509;
2801
2802       fC[12][0] = 0.005;
2803       fC[12][1] = 0.005;
2804       fC[12][2] = 1.0000;
2805       fC[12][3] = 0.3988;
2806       fC[12][4] = 0.5339;
2807
2808       fC[13][0] = 0.005;
2809       fC[13][1] = 0.005;
2810       fC[13][2] = 1.0000;
2811       fC[13][3] = 0.4315;
2812       fC[13][4] = 0.5995;
2813
2814       fC[14][0] = 0.005;
2815       fC[14][1] = 0.005;
2816       fC[14][2] = 1.0000;
2817       fC[14][3] = 0.4548;
2818       fC[14][4] = 0.6612;
2819
2820       fC[15][0] = 0.005;
2821       fC[15][1] = 0.005;
2822       fC[15][2] = 1.0000;
2823       fC[15][3] = 0.4744;
2824       fC[15][4] = 0.7060;
2825
2826       fC[16][0] = 0.005;
2827       fC[16][1] = 0.005;
2828       fC[16][2] = 1.0000;
2829       fC[16][3] = 0.4899;
2830       fC[16][4] = 0.7388;
2831
2832       fC[17][0] = 0.005;
2833       fC[17][1] = 0.005;
2834       fC[17][2] = 1.0000;
2835       fC[17][3] = 0.4411;
2836       fC[17][4] = 0.7293;
2837   }
2838   // 30-40%
2839   else if(centrCur < 40){
2840       fC[0][0] = 0.005;
2841       fC[0][1] = 0.005;
2842       fC[0][2] = 1.0000;
2843       fC[0][3] = 0.0010;
2844       fC[0][4] = 0.0010;
2845
2846       fC[1][0] = 0.005;
2847       fC[1][1] = 0.005;
2848       fC[1][2] = 1.0000;
2849       fC[1][3] = 0.0102;
2850       fC[1][4] = 0.0048;
2851
2852       fC[2][0] = 0.005;
2853       fC[2][1] = 0.005;
2854       fC[2][2] = 1.0000;
2855       fC[2][3] = 0.0306;
2856       fC[2][4] = 0.0079;
2857
2858       fC[3][0] = 0.005;
2859       fC[3][1] = 0.005;
2860       fC[3][2] = 1.0000;
2861       fC[3][3] = 0.0617;
2862       fC[3][4] = 0.0338;
2863
2864       fC[4][0] = 0.005;
2865       fC[4][1] = 0.005;
2866       fC[4][2] = 1.0000;
2867       fC[4][3] = 0.0920;
2868       fC[4][4] = 0.0652;
2869
2870       fC[5][0] = 0.005;
2871       fC[5][1] = 0.005;
2872       fC[5][2] = 1.0000;
2873       fC[5][3] = 0.1211;
2874       fC[5][4] = 0.0955;
2875
2876       fC[6][0] = 0.005;
2877       fC[6][1] = 0.005;
2878       fC[6][2] = 1.0000;
2879       fC[6][3] = 0.1496;
2880       fC[6][4] = 0.1242;
2881
2882       fC[7][0] = 0.005;
2883       fC[7][1] = 0.005;
2884       fC[7][2] = 1.0000;
2885       fC[7][3] = 0.1807;
2886       fC[7][4] = 0.1576;
2887
2888       fC[8][0] = 0.005;
2889       fC[8][1] = 0.005;
2890       fC[8][2] = 1.0000;
2891       fC[8][3] = 0.2195;
2892       fC[8][4] = 0.2097;
2893
2894       fC[9][0] = 0.005;
2895       fC[9][1] = 0.005;
2896       fC[9][2] = 1.0000;
2897       fC[9][3] = 0.2732;
2898       fC[9][4] = 0.2884;
2899
2900       fC[10][0] = 0.005;
2901       fC[10][1] = 0.005;
2902       fC[10][2] = 1.0000;
2903       fC[10][3] = 0.3204;
2904       fC[10][4] = 0.3679;
2905
2906       fC[11][0] = 0.005;
2907       fC[11][1] = 0.005;
2908       fC[11][2] = 1.0000;
2909       fC[11][3] = 0.3564;
2910       fC[11][4] = 0.4449;
2911
2912       fC[12][0] = 0.005;
2913       fC[12][1] = 0.005;
2914       fC[12][2] = 1.0000;
2915       fC[12][3] = 0.3791;
2916       fC[12][4] = 0.5052;
2917
2918       fC[13][0] = 0.005;
2919       fC[13][1] = 0.005;
2920       fC[13][2] = 1.0000;
2921       fC[13][3] = 0.4062;
2922       fC[13][4] = 0.5647;
2923
2924       fC[14][0] = 0.005;
2925       fC[14][1] = 0.005;
2926       fC[14][2] = 1.0000;
2927       fC[14][3] = 0.4234;
2928       fC[14][4] = 0.6203;
2929
2930       fC[15][0] = 0.005;
2931       fC[15][1] = 0.005;
2932       fC[15][2] = 1.0000;
2933       fC[15][3] = 0.4441;
2934       fC[15][4] = 0.6381;
2935
2936       fC[16][0] = 0.005;
2937       fC[16][1] = 0.005;
2938       fC[16][2] = 1.0000;
2939       fC[16][3] = 0.4629;
2940       fC[16][4] = 0.6496;
2941
2942       fC[17][0] = 0.005;
2943       fC[17][1] = 0.005;
2944       fC[17][2] = 1.0000;
2945       fC[17][3] = 0.4293;
2946       fC[17][4] = 0.6491;
2947   }
2948   // 40-50%
2949   else if(centrCur < 50){
2950       fC[0][0] = 0.005;
2951       fC[0][1] = 0.005;
2952       fC[0][2] = 1.0000;
2953       fC[0][3] = 0.0010;
2954       fC[0][4] = 0.0010;
2955
2956       fC[1][0] = 0.005;
2957       fC[1][1] = 0.005;
2958       fC[1][2] = 1.0000;
2959       fC[1][3] = 0.0093;
2960       fC[1][4] = 0.0057;
2961
2962       fC[2][0] = 0.005;
2963       fC[2][1] = 0.005;
2964       fC[2][2] = 1.0000;
2965       fC[2][3] = 0.0319;
2966       fC[2][4] = 0.0075;
2967
2968       fC[3][0] = 0.005;
2969       fC[3][1] = 0.005;
2970       fC[3][2] = 1.0000;
2971       fC[3][3] = 0.0639;
2972       fC[3][4] = 0.0371;
2973
2974       fC[4][0] = 0.005;
2975       fC[4][1] = 0.005;
2976       fC[4][2] = 1.0000;
2977       fC[4][3] = 0.0939;
2978       fC[4][4] = 0.0725;
2979
2980       fC[5][0] = 0.005;
2981       fC[5][1] = 0.005;
2982       fC[5][2] = 1.0000;
2983       fC[5][3] = 0.1224;
2984       fC[5][4] = 0.1045;
2985
2986       fC[6][0] = 0.005;
2987       fC[6][1] = 0.005;
2988       fC[6][2] = 1.0000;
2989       fC[6][3] = 0.1520;
2990       fC[6][4] = 0.1387;
2991
2992       fC[7][0] = 0.005;
2993       fC[7][1] = 0.005;
2994       fC[7][2] = 1.0000;
2995       fC[7][3] = 0.1783;
2996       fC[7][4] = 0.1711;
2997
2998       fC[8][0] = 0.005;
2999       fC[8][1] = 0.005;
3000       fC[8][2] = 1.0000;
3001       fC[8][3] = 0.2202;
3002       fC[8][4] = 0.2269;
3003
3004       fC[9][0] = 0.005;
3005       fC[9][1] = 0.005;
3006       fC[9][2] = 1.0000;
3007       fC[9][3] = 0.2672;
3008       fC[9][4] = 0.2955;
3009
3010       fC[10][0] = 0.005;
3011       fC[10][1] = 0.005;
3012       fC[10][2] = 1.0000;
3013       fC[10][3] = 0.3191;
3014       fC[10][4] = 0.3676;
3015
3016       fC[11][0] = 0.005;
3017       fC[11][1] = 0.005;
3018       fC[11][2] = 1.0000;
3019       fC[11][3] = 0.3434;
3020       fC[11][4] = 0.4321;
3021
3022       fC[12][0] = 0.005;
3023       fC[12][1] = 0.005;
3024       fC[12][2] = 1.0000;
3025       fC[12][3] = 0.3692;
3026       fC[12][4] = 0.4879;
3027
3028       fC[13][0] = 0.005;
3029       fC[13][1] = 0.005;
3030       fC[13][2] = 1.0000;
3031       fC[13][3] = 0.3993;
3032       fC[13][4] = 0.5377;
3033
3034       fC[14][0] = 0.005;
3035       fC[14][1] = 0.005;
3036       fC[14][2] = 1.0000;
3037       fC[14][3] = 0.3818;
3038       fC[14][4] = 0.5547;
3039
3040       fC[15][0] = 0.005;
3041       fC[15][1] = 0.005;
3042       fC[15][2] = 1.0000;
3043       fC[15][3] = 0.4003;
3044       fC[15][4] = 0.5484;
3045
3046       fC[16][0] = 0.005;
3047       fC[16][1] = 0.005;
3048       fC[16][2] = 1.0000;
3049       fC[16][3] = 0.4281;
3050       fC[16][4] = 0.5383;
3051
3052       fC[17][0] = 0.005;
3053       fC[17][1] = 0.005;
3054       fC[17][2] = 1.0000;
3055       fC[17][3] = 0.3960;
3056       fC[17][4] = 0.5374;
3057   }
3058   // 50-60%
3059   else if(centrCur < 60){
3060       fC[0][0] = 0.005;
3061       fC[0][1] = 0.005;
3062       fC[0][2] = 1.0000;
3063       fC[0][3] = 0.0010;
3064       fC[0][4] = 0.0010;
3065
3066       fC[1][0] = 0.005;
3067       fC[1][1] = 0.005;
3068       fC[1][2] = 1.0000;
3069       fC[1][3] = 0.0076;
3070       fC[1][4] = 0.0032;
3071
3072       fC[2][0] = 0.005;
3073       fC[2][1] = 0.005;
3074       fC[2][2] = 1.0000;
3075       fC[2][3] = 0.0329;
3076       fC[2][4] = 0.0085;
3077
3078       fC[3][0] = 0.005;
3079       fC[3][1] = 0.005;
3080       fC[3][2] = 1.0000;
3081       fC[3][3] = 0.0653;
3082       fC[3][4] = 0.0423;
3083
3084       fC[4][0] = 0.005;
3085       fC[4][1] = 0.005;
3086       fC[4][2] = 1.0000;
3087       fC[4][3] = 0.0923;
3088       fC[4][4] = 0.0813;
3089
3090       fC[5][0] = 0.005;
3091       fC[5][1] = 0.005;
3092       fC[5][2] = 1.0000;
3093       fC[5][3] = 0.1219;
3094       fC[5][4] = 0.1161;
3095
3096       fC[6][0] = 0.005;
3097       fC[6][1] = 0.005;
3098       fC[6][2] = 1.0000;
3099       fC[6][3] = 0.1519;
3100       fC[6][4] = 0.1520;
3101
3102       fC[7][0] = 0.005;
3103       fC[7][1] = 0.005;
3104       fC[7][2] = 1.0000;
3105       fC[7][3] = 0.1763;
3106       fC[7][4] = 0.1858;
3107
3108       fC[8][0] = 0.005;
3109       fC[8][1] = 0.005;
3110       fC[8][2] = 1.0000;
3111       fC[8][3] = 0.2178;
3112       fC[8][4] = 0.2385;
3113
3114       fC[9][0] = 0.005;
3115       fC[9][1] = 0.005;
3116       fC[9][2] = 1.0000;
3117       fC[9][3] = 0.2618;
3118       fC[9][4] = 0.3070;
3119
3120       fC[10][0] = 0.005;
3121       fC[10][1] = 0.005;
3122       fC[10][2] = 1.0000;
3123       fC[10][3] = 0.3067;
3124       fC[10][4] = 0.3625;
3125
3126       fC[11][0] = 0.005;
3127       fC[11][1] = 0.005;
3128       fC[11][2] = 1.0000;
3129       fC[11][3] = 0.3336;
3130       fC[11][4] = 0.4188;
3131
3132       fC[12][0] = 0.005;
3133       fC[12][1] = 0.005;
3134       fC[12][2] = 1.0000;
3135       fC[12][3] = 0.3706;
3136       fC[12][4] = 0.4511;
3137
3138       fC[13][0] = 0.005;
3139       fC[13][1] = 0.005;
3140       fC[13][2] = 1.0000;
3141       fC[13][3] = 0.3765;
3142       fC[13][4] = 0.4729;
3143
3144       fC[14][0] = 0.005;
3145       fC[14][1] = 0.005;
3146       fC[14][2] = 1.0000;
3147       fC[14][3] = 0.3942;
3148       fC[14][4] = 0.4855;
3149
3150       fC[15][0] = 0.005;
3151       fC[15][1] = 0.005;
3152       fC[15][2] = 1.0000;
3153       fC[15][3] = 0.4051;
3154       fC[15][4] = 0.4762;
3155
3156       fC[16][0] = 0.005;
3157       fC[16][1] = 0.005;
3158       fC[16][2] = 1.0000;
3159       fC[16][3] = 0.3843;
3160       fC[16][4] = 0.4763;
3161
3162       fC[17][0] = 0.005;
3163       fC[17][1] = 0.005;
3164       fC[17][2] = 1.0000;
3165       fC[17][3] = 0.4237;
3166       fC[17][4] = 0.4773;
3167   }
3168   // 60-70%
3169   else if(centrCur < 70){
3170          fC[0][0] = 0.005;
3171       fC[0][1] = 0.005;
3172       fC[0][2] = 1.0000;
3173       fC[0][3] = 0.0010;
3174       fC[0][4] = 0.0010;
3175
3176       fC[1][0] = 0.005;
3177       fC[1][1] = 0.005;
3178       fC[1][2] = 1.0000;
3179       fC[1][3] = 0.0071;
3180       fC[1][4] = 0.0012;
3181
3182       fC[2][0] = 0.005;
3183       fC[2][1] = 0.005;
3184       fC[2][2] = 1.0000;
3185       fC[2][3] = 0.0336;
3186       fC[2][4] = 0.0097;
3187
3188       fC[3][0] = 0.005;
3189       fC[3][1] = 0.005;
3190       fC[3][2] = 1.0000;
3191       fC[3][3] = 0.0662;
3192       fC[3][4] = 0.0460;
3193
3194       fC[4][0] = 0.005;
3195       fC[4][1] = 0.005;
3196       fC[4][2] = 1.0000;
3197       fC[4][3] = 0.0954;
3198       fC[4][4] = 0.0902;
3199
3200       fC[5][0] = 0.005;
3201       fC[5][1] = 0.005;
3202       fC[5][2] = 1.0000;
3203       fC[5][3] = 0.1181;
3204       fC[5][4] = 0.1306;
3205
3206       fC[6][0] = 0.005;
3207       fC[6][1] = 0.005;
3208       fC[6][2] = 1.0000;
3209       fC[6][3] = 0.1481;
3210       fC[6][4] = 0.1662;
3211
3212       fC[7][0] = 0.005;
3213       fC[7][1] = 0.005;
3214       fC[7][2] = 1.0000;
3215       fC[7][3] = 0.1765;
3216       fC[7][4] = 0.1963;
3217
3218       fC[8][0] = 0.005;
3219       fC[8][1] = 0.005;
3220       fC[8][2] = 1.0000;
3221       fC[8][3] = 0.2155;
3222       fC[8][4] = 0.2433;
3223
3224       fC[9][0] = 0.005;
3225       fC[9][1] = 0.005;
3226       fC[9][2] = 1.0000;
3227       fC[9][3] = 0.2580;
3228       fC[9][4] = 0.3022;
3229
3230       fC[10][0] = 0.005;
3231       fC[10][1] = 0.005;
3232       fC[10][2] = 1.0000;
3233       fC[10][3] = 0.2872;
3234       fC[10][4] = 0.3481;
3235
3236       fC[11][0] = 0.005;
3237       fC[11][1] = 0.005;
3238       fC[11][2] = 1.0000;
3239       fC[11][3] = 0.3170;
3240       fC[11][4] = 0.3847;
3241
3242       fC[12][0] = 0.005;
3243       fC[12][1] = 0.005;
3244       fC[12][2] = 1.0000;
3245       fC[12][3] = 0.3454;
3246       fC[12][4] = 0.4258;
3247
3248       fC[13][0] = 0.005;
3249       fC[13][1] = 0.005;
3250       fC[13][2] = 1.0000;
3251       fC[13][3] = 0.3580;
3252       fC[13][4] = 0.4299;
3253
3254       fC[14][0] = 0.005;
3255       fC[14][1] = 0.005;
3256       fC[14][2] = 1.0000;
3257       fC[14][3] = 0.3903;
3258       fC[14][4] = 0.4326;
3259
3260       fC[15][0] = 0.005;
3261       fC[15][1] = 0.005;
3262       fC[15][2] = 1.0000;
3263       fC[15][3] = 0.3690;
3264       fC[15][4] = 0.4491;
3265
3266       fC[16][0] = 0.005;
3267       fC[16][1] = 0.005;
3268       fC[16][2] = 1.0000;
3269       fC[16][3] = 0.4716;
3270       fC[16][4] = 0.4298;
3271
3272       fC[17][0] = 0.005;
3273       fC[17][1] = 0.005;
3274       fC[17][2] = 1.0000;
3275       fC[17][3] = 0.3875;
3276       fC[17][4] = 0.4083;
3277   }
3278   // 70-80%
3279   else if(centrCur < 80){
3280       fC[0][0] = 0.005;
3281       fC[0][1] = 0.005;
3282       fC[0][2] = 1.0000;
3283       fC[0][3] = 0.0010;
3284       fC[0][4] = 0.0010;
3285
3286       fC[1][0] = 0.005;
3287       fC[1][1] = 0.005;
3288       fC[1][2] = 1.0000;
3289       fC[1][3] = 0.0075;
3290       fC[1][4] = 0.0007;
3291
3292       fC[2][0] = 0.005;
3293       fC[2][1] = 0.005;
3294       fC[2][2] = 1.0000;
3295       fC[2][3] = 0.0313;
3296       fC[2][4] = 0.0124;
3297
3298       fC[3][0] = 0.005;
3299       fC[3][1] = 0.005;
3300       fC[3][2] = 1.0000;
3301       fC[3][3] = 0.0640;
3302       fC[3][4] = 0.0539;
3303
3304       fC[4][0] = 0.005;
3305       fC[4][1] = 0.005;
3306       fC[4][2] = 1.0000;
3307       fC[4][3] = 0.0923;
3308       fC[4][4] = 0.0992;
3309
3310       fC[5][0] = 0.005;
3311       fC[5][1] = 0.005;
3312       fC[5][2] = 1.0000;
3313       fC[5][3] = 0.1202;
3314       fC[5][4] = 0.1417;
3315
3316       fC[6][0] = 0.005;
3317       fC[6][1] = 0.005;
3318       fC[6][2] = 1.0000;
3319       fC[6][3] = 0.1413;
3320       fC[6][4] = 0.1729;
3321
3322       fC[7][0] = 0.005;
3323       fC[7][1] = 0.005;
3324       fC[7][2] = 1.0000;
3325       fC[7][3] = 0.1705;
3326       fC[7][4] = 0.1999;
3327
3328       fC[8][0] = 0.005;
3329       fC[8][1] = 0.005;
3330       fC[8][2] = 1.0000;
3331       fC[8][3] = 0.2103;
3332       fC[8][4] = 0.2472;
3333
3334       fC[9][0] = 0.005;
3335       fC[9][1] = 0.005;
3336       fC[9][2] = 1.0000;
3337       fC[9][3] = 0.2373;
3338       fC[9][4] = 0.2916;
3339
3340       fC[10][0] = 0.005;
3341       fC[10][1] = 0.005;
3342       fC[10][2] = 1.0000;
3343       fC[10][3] = 0.2824;
3344       fC[10][4] = 0.3323;
3345
3346       fC[11][0] = 0.005;
3347       fC[11][1] = 0.005;
3348       fC[11][2] = 1.0000;
3349       fC[11][3] = 0.3046;
3350       fC[11][4] = 0.3576;
3351
3352       fC[12][0] = 0.005;
3353       fC[12][1] = 0.005;
3354       fC[12][2] = 1.0000;
3355       fC[12][3] = 0.3585;
3356       fC[12][4] = 0.4003;
3357
3358       fC[13][0] = 0.005;
3359       fC[13][1] = 0.005;
3360       fC[13][2] = 1.0000;
3361       fC[13][3] = 0.3461;
3362       fC[13][4] = 0.3982;
3363
3364       fC[14][0] = 0.005;
3365       fC[14][1] = 0.005;
3366       fC[14][2] = 1.0000;
3367       fC[14][3] = 0.3362;
3368       fC[14][4] = 0.3776;
3369
3370       fC[15][0] = 0.005;
3371       fC[15][1] = 0.005;
3372       fC[15][2] = 1.0000;
3373       fC[15][3] = 0.3071;
3374       fC[15][4] = 0.3500;
3375
3376       fC[16][0] = 0.005;
3377       fC[16][1] = 0.005;
3378       fC[16][2] = 1.0000;
3379       fC[16][3] = 0.2914;
3380       fC[16][4] = 0.3937;
3381
3382       fC[17][0] = 0.005;
3383       fC[17][1] = 0.005;
3384       fC[17][2] = 1.0000;
3385       fC[17][3] = 0.3727;
3386       fC[17][4] = 0.3877;
3387   }
3388   // 80-100%
3389   else{
3390       fC[0][0] = 0.005;
3391       fC[0][1] = 0.005;
3392       fC[0][2] = 1.0000;
3393       fC[0][3] = 0.0010;
3394       fC[0][4] = 0.0010;
3395
3396       fC[1][0] = 0.005;
3397       fC[1][1] = 0.005;
3398       fC[1][2] = 1.0000;
3399       fC[1][3] = 0.0060;
3400       fC[1][4] = 0.0035;
3401
3402       fC[2][0] = 0.005;
3403       fC[2][1] = 0.005;
3404       fC[2][2] = 1.0000;
3405       fC[2][3] = 0.0323;
3406       fC[2][4] = 0.0113;
3407
3408       fC[3][0] = 0.005;
3409       fC[3][1] = 0.005;
3410       fC[3][2] = 1.0000;
3411       fC[3][3] = 0.0609;
3412       fC[3][4] = 0.0653;
3413
3414       fC[4][0] = 0.005;
3415       fC[4][1] = 0.005;
3416       fC[4][2] = 1.0000;
3417       fC[4][3] = 0.0922;
3418       fC[4][4] = 0.1076;
3419
3420       fC[5][0] = 0.005;
3421       fC[5][1] = 0.005;
3422       fC[5][2] = 1.0000;
3423       fC[5][3] = 0.1096;
3424       fC[5][4] = 0.1328;
3425
3426       fC[6][0] = 0.005;
3427       fC[6][1] = 0.005;
3428       fC[6][2] = 1.0000;
3429       fC[6][3] = 0.1495;
3430       fC[6][4] = 0.1779;
3431
3432       fC[7][0] = 0.005;
3433       fC[7][1] = 0.005;
3434       fC[7][2] = 1.0000;
3435       fC[7][3] = 0.1519;
3436       fC[7][4] = 0.1989;
3437
3438       fC[8][0] = 0.005;
3439       fC[8][1] = 0.005;
3440       fC[8][2] = 1.0000;
3441       fC[8][3] = 0.1817;
3442       fC[8][4] = 0.2472;
3443
3444       fC[9][0] = 0.005;
3445       fC[9][1] = 0.005;
3446       fC[9][2] = 1.0000;
3447       fC[9][3] = 0.2429;
3448       fC[9][4] = 0.2684;
3449
3450       fC[10][0] = 0.005;
3451       fC[10][1] = 0.005;
3452       fC[10][2] = 1.0000;
3453       fC[10][3] = 0.2760;
3454       fC[10][4] = 0.3098;
3455
3456       fC[11][0] = 0.005;
3457       fC[11][1] = 0.005;
3458       fC[11][2] = 1.0000;
3459       fC[11][3] = 0.2673;
3460       fC[11][4] = 0.3198;
3461
3462       fC[12][0] = 0.005;
3463       fC[12][1] = 0.005;
3464       fC[12][2] = 1.0000;
3465       fC[12][3] = 0.3165;
3466       fC[12][4] = 0.3564;
3467
3468       fC[13][0] = 0.005;
3469       fC[13][1] = 0.005;
3470       fC[13][2] = 1.0000;
3471       fC[13][3] = 0.3526;
3472       fC[13][4] = 0.3011;
3473
3474       fC[14][0] = 0.005;
3475       fC[14][1] = 0.005;
3476       fC[14][2] = 1.0000;
3477       fC[14][3] = 0.3788;
3478       fC[14][4] = 0.3011;
3479
3480       fC[15][0] = 0.005;
3481       fC[15][1] = 0.005;
3482       fC[15][2] = 1.0000;
3483       fC[15][3] = 0.3788;
3484       fC[15][4] = 0.3011;
3485
3486       fC[16][0] = 0.005;
3487       fC[16][1] = 0.005;
3488       fC[16][2] = 1.0000;
3489       fC[16][3] = 0.3788;
3490       fC[16][4] = 0.3011;
3491
3492       fC[17][0] = 0.005;
3493       fC[17][1] = 0.005;
3494       fC[17][2] = 1.0000;
3495       fC[17][3] = 0.3788;
3496       fC[17][4] = 0.3011;
3497   }
3498
3499   for(Int_t i=18;i<fgkPIDptBin;i++){
3500     fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
3501     fC[i][0] = fC[17][0];
3502     fC[i][1] = fC[17][1];
3503     fC[i][2] = fC[17][2];
3504     fC[i][3] = fC[17][3];
3505     fC[i][4] = fC[17][4];
3506   }
3507 }
3508
3509 //---------------------------------------------------------------//
3510 Bool_t AliFlowTrackCuts::TPCTOFagree(const AliESDtrack *track)
3511 {
3512   //check pid agreement between TPC and TOF
3513   Bool_t status = kFALSE;
3514
3515   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
3516
3517
3518   Double_t exptimes[5];
3519   track->GetIntegratedTimes(exptimes);
3520
3521   Float_t dedx = track->GetTPCsignal();
3522
3523   Float_t p = track->P();
3524   Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
3525   Float_t tl = track->GetIntegratedLength();
3526
3527   Float_t betagammares =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
3528
3529   Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
3530
3531 //  printf("betagamma1 = %f\n",betagamma1);
3532
3533   if(betagamma1 < 0.1) betagamma1 = 0.1;
3534
3535   if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
3536   else betagamma1 = 100;
3537
3538   Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
3539 //  printf("betagamma2 = %f\n",betagamma2);
3540
3541   if(betagamma2 < 0.1) betagamma2 = 0.1;
3542
3543   if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
3544   else betagamma2 = 100;
3545
3546
3547   Double_t ptpc[3];
3548   track->GetInnerPxPyPz(ptpc);
3549   Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
3550
3551   for(Int_t i=0;i < 5;i++){
3552     Float_t resolutionTOF =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
3553     if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
3554       Float_t dedxExp = 0;
3555       if(i==0) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
3556       else if(i==1) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
3557       else if(i==2) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
3558       else if(i==3) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
3559       else if(i==4) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
3560
3561       Float_t resolutionTPC = 2;
3562       if(i==0) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
3563       else if(i==1) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
3564       else if(i==2) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
3565       else if(i==3) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
3566       else if(i==4) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
3567
3568       if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
3569         status = kTRUE;
3570       }
3571     }
3572   }
3573
3574   Float_t bb1 =  fESDpid.GetTPCResponse().Bethe(betagamma1);
3575   Float_t bb2 =  fESDpid.GetTPCResponse().Bethe(betagamma2);
3576   Float_t bbM =  fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
3577
3578
3579   //  status = kFALSE;
3580   // for nuclei
3581   Float_t resolutionTOFpr =   fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
3582   Float_t resolutionTPCpr =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
3583   if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3584      status = kTRUE;
3585   }
3586   else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3587      status = kTRUE;
3588   }
3589   else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3590      status = kTRUE;
3591   }
3592
3593   return status;
3594 }
3595 // end part added by F. Noferini
3596 //-----------------------------------------------------------------------
3597
3598 //-----------------------------------------------------------------------
3599 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
3600 {
3601   //get the name of the particle id source
3602   switch (s)
3603   {
3604     case kTPCdedx:
3605       return "TPCdedx";
3606     case kTPCbayesian:
3607       return "TPCbayesian";
3608     case kTOFbeta:
3609       return "TOFbeta";
3610     case kTPCpid:
3611       return "TPCpid";
3612     case kTOFpid:
3613       return "TOFpid";
3614     case kTOFbayesian:
3615       return "TOFbayesianPID";
3616     case kTOFbetaSimple:
3617       return "TOFbetaSimple";
3618     case kTPCNuclei:
3619       return "TPCnuclei";
3620     default:
3621       return "NOPID";
3622   }
3623 }
3624
3625 //-----------------------------------------------------------------------
3626 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
3627 {
3628   //return the name of the selected parameter type
3629   switch (type)
3630   {
3631     case kMC:
3632       return "MC";
3633     case kGlobal:
3634       return "Global";
3635     case kTPCstandalone:
3636       return "TPCstandalone";
3637     case kSPDtracklet:
3638       return "SPDtracklets";
3639     case kPMD:
3640       return "PMD";
3641     case kV0:
3642       return "V0";
3643     default:
3644       return "unknown";
3645   }
3646 }
3647
3648 //-----------------------------------------------------------------------
3649 Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
3650 {
3651   //check PMD specific cuts
3652   //clean up from last iteration, and init label
3653   Int_t   det   = track->GetDetector();
3654   //Int_t   smn   = track->GetSmn();
3655   Float_t clsX  = track->GetClusterX();
3656   Float_t clsY  = track->GetClusterY();
3657   Float_t clsZ  = track->GetClusterZ();
3658   Float_t ncell = track->GetClusterCells();
3659   Float_t adc   = track->GetClusterADC();
3660
3661   fTrack = NULL;
3662   fMCparticle=NULL;
3663   fTrackLabel=-996;
3664
3665   fTrackEta = GetPmdEta(clsX,clsY,clsZ);
3666   fTrackPhi = GetPmdPhi(clsX,clsY);
3667   fTrackWeight = 1.0;
3668
3669   Bool_t pass=kTRUE;
3670   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
3671   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
3672   if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
3673   if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
3674   if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
3675
3676   return pass;
3677 }
3678
3679 //-----------------------------------------------------------------------
3680 Bool_t AliFlowTrackCuts::PassesV0cuts(Int_t id)
3681 {
3682   //check V0 cuts
3683   if (id<0) return kFALSE;
3684
3685   //clean up from last iter
3686   fTrack = NULL;
3687   fMCparticle=NULL;
3688   fTrackLabel=-995;
3689
3690   fTrackPhi = TMath::PiOver4()*(0.5+id%8);
3691
3692   if(id<32)
3693     fTrackEta = -3.45+0.5*(id/8);
3694   else
3695     fTrackEta = +4.8-0.6*((id/8)-4);
3696
3697   fTrackWeight = fEvent->GetVZEROEqMultiplicity(id);
3698
3699   if (fLinearizeVZEROresponse)
3700   {
3701     //this is only needed in pass1 of LHC10h
3702     Float_t multV0[fgkNumberOfV0tracks];
3703     Float_t dummy=0.0;
3704     AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multV0);
3705     fTrackWeight = multV0[id];
3706   }
3707
3708   Bool_t pass=kTRUE;
3709   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
3710   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
3711
3712   return pass;
3713 }
3714
3715 //----------------------------------------------------------------------------//
3716 Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
3717 {
3718   //get PMD "track" eta coordinate
3719   Float_t rpxpy, theta, eta;
3720   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
3721   theta  = TMath::ATan2(rpxpy,zPos);
3722   eta    = -TMath::Log(TMath::Tan(0.5*theta));
3723   return eta;
3724 }
3725
3726 //--------------------------------------------------------------------------//
3727 Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
3728 {
3729   //Get PMD "track" phi coordinate
3730   Float_t pybypx, phi = 0., phi1;
3731   if(xPos==0)
3732     {
3733       if(yPos>0) phi = 90.;
3734       if(yPos<0) phi = 270.;
3735     }
3736   if(xPos != 0)
3737     {
3738       pybypx = yPos/xPos;
3739       if(pybypx < 0) pybypx = - pybypx;
3740       phi1 = TMath::ATan(pybypx)*180./3.14159;
3741
3742       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
3743       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
3744       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
3745       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
3746
3747     }
3748   phi = phi*3.14159/180.;
3749   return   phi;
3750 }
3751
3752 //---------------------------------------------------------------//
3753 void AliFlowTrackCuts::Browse(TBrowser* b)
3754 {
3755   //some browsing capabilities
3756   if (fQA) b->Add(fQA);
3757 }
3758
3759 //---------------------------------------------------------------//
3760 Long64_t AliFlowTrackCuts::Merge(TCollection* list)
3761 {
3762   //merge
3763   if (!fQA || !list) return 0;
3764   if (list->IsEmpty()) return 0;
3765   AliFlowTrackCuts* obj=NULL;
3766   TList tmplist;
3767   TIter next(list);
3768   while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
3769   {
3770     if (obj==this) continue;
3771     tmplist.Add(obj->GetQA());
3772   }
3773   return fQA->Merge(&tmplist);
3774 }
3775
3776
3777