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