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