PWG/FLOW/Tasks/AliFlowTrackCuts.cxx

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