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