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, Int_t passN)
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   cuts->SetMuonPassNumber(passN);
1214   return cuts;
1215 }
1216
1217 //-----------------------------------------------------------------------
1218 Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
1219 {
1220   //fill a flow track from tracklet,vzero,pmd,...
1221   TParticle *tmpTParticle=NULL;
1222   AliMCParticle* tmpAliMCParticle=NULL;
1223   switch (fParamMix)
1224   {
1225     case kPure:
1226       flowtrack->SetPhi(fTrackPhi);
1227       flowtrack->SetEta(fTrackEta);
1228       flowtrack->SetWeight(fTrackWeight);
1229       break;
1230     case kTrackWithMCkine:
1231       if (!fMCparticle) return kFALSE;
1232       flowtrack->SetPhi( fMCparticle->Phi() );
1233       flowtrack->SetEta( fMCparticle->Eta() );
1234       flowtrack->SetPt( fMCparticle->Pt() );
1235       break;
1236     case kTrackWithMCpt:
1237       if (!fMCparticle) return kFALSE;
1238       flowtrack->SetPhi(fTrackPhi);
1239       flowtrack->SetEta(fTrackEta);
1240       flowtrack->SetWeight(fTrackWeight);
1241       flowtrack->SetPt(fMCparticle->Pt());
1242       break;
1243     case kTrackWithPtFromFirstMother:
1244       if (!fMCparticle) return kFALSE;
1245       flowtrack->SetPhi(fTrackPhi);
1246       flowtrack->SetEta(fTrackEta);
1247       flowtrack->SetWeight(fTrackWeight);
1248       tmpTParticle = fMCparticle->Particle();
1249       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1250       flowtrack->SetPt(tmpAliMCParticle->Pt());
1251       break;
1252     default:
1253       flowtrack->SetPhi(fTrackPhi);
1254       flowtrack->SetEta(fTrackEta);
1255       flowtrack->SetWeight(fTrackWeight);
1256       break;
1257   }
1258   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1259   return kTRUE;
1260 }
1261
1262 //-----------------------------------------------------------------------
1263 Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
1264 {
1265   //fill flow track from AliVParticle (ESD,AOD,MC)
1266   if (!fTrack) return kFALSE;
1267   TParticle *tmpTParticle=NULL;
1268   AliMCParticle* tmpAliMCParticle=NULL;
1269   AliExternalTrackParam* externalParams=NULL;
1270   AliESDtrack* esdtrack=NULL;
1271   switch(fParamMix)
1272   {
1273     case kPure:
1274       flowtrack->Set(fTrack);
1275       break;
1276     case kTrackWithMCkine:
1277       flowtrack->Set(fMCparticle);
1278       break;
1279     case kTrackWithMCPID:
1280       flowtrack->Set(fTrack);
1281       //flowtrack->setPID(...) from mc, when implemented
1282       break;
1283     case kTrackWithMCpt:
1284       if (!fMCparticle) return kFALSE;
1285       flowtrack->Set(fTrack);
1286       flowtrack->SetPt(fMCparticle->Pt());
1287       break;
1288     case kTrackWithPtFromFirstMother:
1289       if (!fMCparticle) return kFALSE;
1290       flowtrack->Set(fTrack);
1291       tmpTParticle = fMCparticle->Particle();
1292       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1293       flowtrack->SetPt(tmpAliMCParticle->Pt());
1294       break;
1295     case kTrackWithTPCInnerParams:
1296       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1297       if (!esdtrack) return kFALSE;
1298       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1299       if (!externalParams) return kFALSE;
1300       flowtrack->Set(externalParams);
1301       break;
1302     default:
1303       flowtrack->Set(fTrack);
1304       break;
1305   }
1306   if (fParamType==kMC)
1307   {
1308     flowtrack->SetSource(AliFlowTrack::kFromMC);
1309     flowtrack->SetID(fTrackLabel);
1310   }
1311   else if (dynamic_cast<AliESDtrack*>(fTrack))
1312   {
1313     flowtrack->SetSource(AliFlowTrack::kFromESD);
1314     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1315   }
1316   else if (dynamic_cast<AliESDMuonTrack*>(fTrack))                                  // XZhang 20120604
1317   {                                                                                 // XZhang 20120604
1318     flowtrack->SetSource(AliFlowTrack::kFromMUON);                                  // XZhang 20120604
1319     flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID());  // XZhang 20120604
1320   }                                                                                 // XZhang 20120604
1321   else if (dynamic_cast<AliAODTrack*>(fTrack))
1322   {
1323     if (fParamType==kMUON)                            // XZhang 20120604
1324       flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1325     else                                              // XZhang 20120604
1326       flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1327     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1328   }
1329   else if (dynamic_cast<AliMCParticle*>(fTrack))
1330   {
1331     flowtrack->SetSource(AliFlowTrack::kFromMC);
1332     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1333   }
1334   return kTRUE;
1335 }
1336
1337 //-----------------------------------------------------------------------
1338 Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
1339 {
1340   //fill a flow track constructed from whatever we applied cuts on
1341   //return true on success
1342   switch (fParamType)
1343   {
1344     case kSPDtracklet:
1345       return FillFlowTrackGeneric(track);
1346     case kPMD:
1347       return FillFlowTrackGeneric(track);
1348     case kV0:
1349       return FillFlowTrackGeneric(track);
1350     default:
1351       return FillFlowTrackVParticle(track);
1352   }
1353 }
1354
1355 //-----------------------------------------------------------------------
1356 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
1357 {
1358   //make a flow track from tracklet
1359   AliFlowTrack* flowtrack=NULL;
1360   TParticle *tmpTParticle=NULL;
1361   AliMCParticle* tmpAliMCParticle=NULL;
1362   switch (fParamMix)
1363   {
1364     case kPure:
1365       flowtrack = new AliFlowTrack();
1366       flowtrack->SetPhi(fTrackPhi);
1367       flowtrack->SetEta(fTrackEta);
1368       flowtrack->SetWeight(fTrackWeight);
1369       break;
1370     case kTrackWithMCkine:
1371       if (!fMCparticle) return NULL;
1372       flowtrack = new AliFlowTrack();
1373       flowtrack->SetPhi( fMCparticle->Phi() );
1374       flowtrack->SetEta( fMCparticle->Eta() );
1375       flowtrack->SetPt( fMCparticle->Pt() );
1376       break;
1377     case kTrackWithMCpt:
1378       if (!fMCparticle) return NULL;
1379       flowtrack = new AliFlowTrack();
1380       flowtrack->SetPhi(fTrackPhi);
1381       flowtrack->SetEta(fTrackEta);
1382       flowtrack->SetWeight(fTrackWeight);
1383       flowtrack->SetPt(fMCparticle->Pt());
1384       break;
1385     case kTrackWithPtFromFirstMother:
1386       if (!fMCparticle) return NULL;
1387       flowtrack = new AliFlowTrack();
1388       flowtrack->SetPhi(fTrackPhi);
1389       flowtrack->SetEta(fTrackEta);
1390       flowtrack->SetWeight(fTrackWeight);
1391       tmpTParticle = fMCparticle->Particle();
1392       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1393       flowtrack->SetPt(tmpAliMCParticle->Pt());
1394       break;
1395     default:
1396       flowtrack = new AliFlowTrack();
1397       flowtrack->SetPhi(fTrackPhi);
1398       flowtrack->SetEta(fTrackEta);
1399       flowtrack->SetWeight(fTrackWeight);
1400       break;
1401   }
1402   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1403   return flowtrack;
1404 }
1405
1406 //-----------------------------------------------------------------------
1407 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
1408 {
1409   //make flow track from AliVParticle (ESD,AOD,MC)
1410   if (!fTrack) return NULL;
1411   AliFlowTrack* flowtrack=NULL;
1412   TParticle *tmpTParticle=NULL;
1413   AliMCParticle* tmpAliMCParticle=NULL;
1414   AliExternalTrackParam* externalParams=NULL;
1415   AliESDtrack* esdtrack=NULL;
1416   switch(fParamMix)
1417   {
1418     case kPure:
1419       flowtrack = new AliFlowTrack(fTrack);
1420       break;
1421     case kTrackWithMCkine:
1422       flowtrack = new AliFlowTrack(fMCparticle);
1423       break;
1424     case kTrackWithMCPID:
1425       flowtrack = new AliFlowTrack(fTrack);
1426       //flowtrack->setPID(...) from mc, when implemented
1427       break;
1428     case kTrackWithMCpt:
1429       if (!fMCparticle) return NULL;
1430       flowtrack = new AliFlowTrack(fTrack);
1431       flowtrack->SetPt(fMCparticle->Pt());
1432       break;
1433     case kTrackWithPtFromFirstMother:
1434       if (!fMCparticle) return NULL;
1435       flowtrack = new AliFlowTrack(fTrack);
1436       tmpTParticle = fMCparticle->Particle();
1437       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1438       flowtrack->SetPt(tmpAliMCParticle->Pt());
1439       break;
1440     case kTrackWithTPCInnerParams:
1441       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1442       if (!esdtrack) return NULL;
1443       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1444       if (!externalParams) return NULL;
1445       flowtrack = new AliFlowTrack(externalParams);
1446       break;
1447     default:
1448       flowtrack = new AliFlowTrack(fTrack);
1449       break;
1450   }
1451   if (fParamType==kMC)
1452   {
1453     flowtrack->SetSource(AliFlowTrack::kFromMC);
1454     flowtrack->SetID(fTrackLabel);
1455   }
1456   else if (dynamic_cast<AliESDtrack*>(fTrack))
1457   {
1458     flowtrack->SetSource(AliFlowTrack::kFromESD);
1459     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1460   }
1461   else if (dynamic_cast<AliAODTrack*>(fTrack))
1462   {
1463     flowtrack->SetSource(AliFlowTrack::kFromAOD);
1464     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1465   }
1466   else if (dynamic_cast<AliMCParticle*>(fTrack))
1467   {
1468     flowtrack->SetSource(AliFlowTrack::kFromMC);
1469     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1470   }
1471   return flowtrack;
1472 }
1473
1474 //-----------------------------------------------------------------------
1475 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
1476 {
1477   //make a flow track from PMD track
1478   AliFlowTrack* flowtrack=NULL;
1479   TParticle *tmpTParticle=NULL;
1480   AliMCParticle* tmpAliMCParticle=NULL;
1481   switch (fParamMix)
1482   {
1483     case kPure:
1484       flowtrack = new AliFlowTrack();
1485       flowtrack->SetPhi(fTrackPhi);
1486       flowtrack->SetEta(fTrackEta);
1487       flowtrack->SetWeight(fTrackWeight);
1488       break;
1489     case kTrackWithMCkine:
1490       if (!fMCparticle) return NULL;
1491       flowtrack = new AliFlowTrack();
1492       flowtrack->SetPhi( fMCparticle->Phi() );
1493       flowtrack->SetEta( fMCparticle->Eta() );
1494       flowtrack->SetWeight(fTrackWeight);
1495       flowtrack->SetPt( fMCparticle->Pt() );
1496       break;
1497     case kTrackWithMCpt:
1498       if (!fMCparticle) return NULL;
1499       flowtrack = new AliFlowTrack();
1500       flowtrack->SetPhi(fTrackPhi);
1501       flowtrack->SetEta(fTrackEta);
1502       flowtrack->SetWeight(fTrackWeight);
1503       flowtrack->SetPt(fMCparticle->Pt());
1504       break;
1505     case kTrackWithPtFromFirstMother:
1506       if (!fMCparticle) return NULL;
1507       flowtrack = new AliFlowTrack();
1508       flowtrack->SetPhi(fTrackPhi);
1509       flowtrack->SetEta(fTrackEta);
1510       flowtrack->SetWeight(fTrackWeight);
1511       tmpTParticle = fMCparticle->Particle();
1512       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1513       flowtrack->SetPt(tmpAliMCParticle->Pt());
1514       break;
1515     default:
1516       flowtrack = new AliFlowTrack();
1517       flowtrack->SetPhi(fTrackPhi);
1518       flowtrack->SetEta(fTrackEta);
1519       flowtrack->SetWeight(fTrackWeight);
1520       break;
1521   }
1522
1523   flowtrack->SetSource(AliFlowTrack::kFromPMD);
1524   return flowtrack;
1525 }
1526
1527 //-----------------------------------------------------------------------
1528 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackV0() const
1529 {
1530   //make a flow track from V0
1531   AliFlowTrack* flowtrack=NULL;
1532   TParticle *tmpTParticle=NULL;
1533   AliMCParticle* tmpAliMCParticle=NULL;
1534   switch (fParamMix)
1535   {
1536     case kPure:
1537       flowtrack = new AliFlowTrack();
1538       flowtrack->SetPhi(fTrackPhi);
1539       flowtrack->SetEta(fTrackEta);
1540       flowtrack->SetWeight(fTrackWeight);
1541       break;
1542     case kTrackWithMCkine:
1543       if (!fMCparticle) return NULL;
1544       flowtrack = new AliFlowTrack();
1545       flowtrack->SetPhi( fMCparticle->Phi() );
1546       flowtrack->SetEta( fMCparticle->Eta() );
1547       flowtrack->SetWeight(fTrackWeight);
1548       flowtrack->SetPt( fMCparticle->Pt() );
1549       break;
1550     case kTrackWithMCpt:
1551       if (!fMCparticle) return NULL;
1552       flowtrack = new AliFlowTrack();
1553       flowtrack->SetPhi(fTrackPhi);
1554       flowtrack->SetEta(fTrackEta);
1555       flowtrack->SetWeight(fTrackWeight);
1556       flowtrack->SetPt(fMCparticle->Pt());
1557       break;
1558     case kTrackWithPtFromFirstMother:
1559       if (!fMCparticle) return NULL;
1560       flowtrack = new AliFlowTrack();
1561       flowtrack->SetPhi(fTrackPhi);
1562       flowtrack->SetEta(fTrackEta);
1563       flowtrack->SetWeight(fTrackWeight);
1564       tmpTParticle = fMCparticle->Particle();
1565       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1566       flowtrack->SetPt(tmpAliMCParticle->Pt());
1567       break;
1568     default:
1569       flowtrack = new AliFlowTrack();
1570       flowtrack->SetPhi(fTrackPhi);
1571       flowtrack->SetEta(fTrackEta);
1572       flowtrack->SetWeight(fTrackWeight);
1573       break;
1574   }
1575
1576   flowtrack->SetSource(AliFlowTrack::kFromV0);
1577   return flowtrack;
1578 }
1579
1580 //-----------------------------------------------------------------------
1581 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
1582 {
1583   //get a flow track constructed from whatever we applied cuts on
1584   //caller is resposible for deletion
1585   //if construction fails return NULL
1586   //TODO: for tracklets, PMD and V0 we probably need just one method,
1587   //something like MakeFlowTrackGeneric(), wait with this until
1588   //requirements quirks are known.
1589   switch (fParamType)
1590   {
1591     case kSPDtracklet:
1592       return MakeFlowTrackSPDtracklet();
1593     case kPMD:
1594       return MakeFlowTrackPMDtrack();
1595     case kV0:
1596       return MakeFlowTrackV0();
1597     default:
1598       return MakeFlowTrackVParticle();
1599   }
1600 }
1601
1602 //-----------------------------------------------------------------------
1603 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
1604 {
1605   //check if current particle is a physical primary
1606   if (!fMCevent) return kFALSE;
1607   if (fTrackLabel<0) return kFALSE;
1608   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
1609 }
1610
1611 //-----------------------------------------------------------------------
1612 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
1613 {
1614   //check if current particle is a physical primary
1615   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
1616   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
1617   if (!track) return kFALSE;
1618   TParticle* particle = track->Particle();
1619   Bool_t transported = particle->TestBit(kTransportBit);
1620   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
1621         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
1622   return (physprim && (transported || !requiretransported));
1623 }
1624
1625 //-----------------------------------------------------------------------
1626 void AliFlowTrackCuts::DefineHistograms()
1627 {
1628   //define qa histograms
1629   if (fQA) return;
1630
1631   const Int_t kNbinsP=200;
1632   Double_t binsP[kNbinsP+1];
1633   binsP[0]=0.0;
1634   for(int i=1; i<kNbinsP+1; i++)
1635   {
1636     //if(binsP[i-1]+0.05<1.01)
1637     //  binsP[i]=binsP[i-1]+0.05;
1638     //else
1639     binsP[i]=binsP[i-1]+0.05;
1640   }
1641
1642   const Int_t nBinsDCA=1000;
1643   Double_t binsDCA[nBinsDCA+1];
1644   for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
1645   //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
1646
1647   Bool_t adddirstatus = TH1::AddDirectoryStatus();
1648   TH1::AddDirectory(kFALSE);
1649   fQA=new TList(); fQA->SetOwner();
1650   fQA->SetName(Form("%s QA",GetName()));
1651   TList* before = new TList(); before->SetOwner();
1652   before->SetName("before");
1653   TList* after = new TList(); after->SetOwner();
1654   after->SetName("after");
1655   fQA->Add(before);
1656   fQA->Add(after);
1657   before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
1658   after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
1659   before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
1660   after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
1661   before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
1662   after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
1663   //primary
1664   TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
1665   TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
1666   TAxis* axis = NULL;
1667   axis = hb->GetYaxis();
1668   axis->SetBinLabel(1,"secondary");
1669   axis->SetBinLabel(2,"primary");
1670   axis = ha->GetYaxis();
1671   axis->SetBinLabel(1,"secondary");
1672   axis->SetBinLabel(2,"primary");
1673   before->Add(hb); //3
1674   after->Add(ha); //3
1675   //production process
1676   hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
1677                       -0.5, kMaxMCProcess-0.5);
1678   ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
1679                       -0.5, kMaxMCProcess-0.5);
1680   axis = hb->GetYaxis();
1681   for (Int_t i=0; i<kMaxMCProcess; i++)
1682   {
1683     axis->SetBinLabel(i+1,TMCProcessName[i]);
1684   }
1685   axis = ha->GetYaxis();
1686   for (Int_t i=0; i<kMaxMCProcess; i++)
1687   {
1688     axis->SetBinLabel(i+1,TMCProcessName[i]);
1689   }
1690   before->Add(hb); //4
1691   after->Add(ha); //4
1692   //DCA
1693   before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
1694   after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
1695   before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
1696   after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
1697   //first mother
1698   hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
1699   ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
1700   hb->GetYaxis()->SetBinLabel(1, "#gamma");
1701   ha->GetYaxis()->SetBinLabel(1, "#gamma");
1702   hb->GetYaxis()->SetBinLabel(2, "e^{+}");
1703   ha->GetYaxis()->SetBinLabel(2, "e^{+}");
1704   hb->GetYaxis()->SetBinLabel(3, "e^{-}");
1705   ha->GetYaxis()->SetBinLabel(3, "e^{-}");
1706   hb->GetYaxis()->SetBinLabel(4, "#nu");
1707   ha->GetYaxis()->SetBinLabel(4, "#nu");
1708   hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
1709   ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
1710   hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
1711   ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
1712   hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
1713   ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
1714   hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
1715   ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
1716   hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
1717   ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
1718   hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
1719   ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
1720   hb->GetYaxis()->SetBinLabel(11, "K^{+}");
1721   ha->GetYaxis()->SetBinLabel(11, "K^{+}");
1722   hb->GetYaxis()->SetBinLabel(12, "K^{-}");
1723   ha->GetYaxis()->SetBinLabel(12, "K^{-}");
1724   hb->GetYaxis()->SetBinLabel( 13, "n");
1725   ha->GetYaxis()->SetBinLabel( 13, "n");
1726   hb->GetYaxis()->SetBinLabel( 14, "p");
1727   ha->GetYaxis()->SetBinLabel( 14, "p");
1728   hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
1729   ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
1730   hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
1731   ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
1732   hb->GetYaxis()->SetBinLabel(17, "#eta");
1733   ha->GetYaxis()->SetBinLabel(17, "#eta");
1734   hb->GetYaxis()->SetBinLabel(18, "#Lambda");
1735   ha->GetYaxis()->SetBinLabel(18, "#Lambda");
1736   hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
1737   ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
1738   hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
1739   ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
1740   hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
1741   ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
1742   hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
1743   ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
1744   hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
1745   ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
1746   hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
1747   ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
1748   hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
1749   ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
1750   hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
1751   ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
1752   hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
1753   ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
1754   hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
1755   ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
1756   hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
1757   ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
1758   hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
1759   ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
1760   hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
1761   ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
1762   hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
1763   ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
1764   hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
1765   ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
1766   hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
1767   ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
1768   hb->GetYaxis()->SetBinLabel(35, "D^{+}");
1769   ha->GetYaxis()->SetBinLabel(35, "D^{+}");
1770   hb->GetYaxis()->SetBinLabel(36, "D^{-}");
1771   ha->GetYaxis()->SetBinLabel(36, "D^{-}");
1772   hb->GetYaxis()->SetBinLabel(37, "D^{0}");
1773   ha->GetYaxis()->SetBinLabel(37, "D^{0}");
1774   hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
1775   ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
1776   hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
1777   ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
1778   hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
1779   ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
1780   hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
1781   ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
1782   hb->GetYaxis()->SetBinLabel(42, "W^{+}");
1783   ha->GetYaxis()->SetBinLabel(42, "W^{+}");
1784   hb->GetYaxis()->SetBinLabel(43, "W^{-}");
1785   ha->GetYaxis()->SetBinLabel(43, "W^{-}");
1786   hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
1787   ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
1788   before->Add(hb);//7
1789   after->Add(ha);//7
1790
1791   before->Add(new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
1792   after->Add( new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
1793
1794   before->Add(new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
1795   after->Add( new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
1796
1797   before->Add(new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
1798   after->Add( new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
1799
1800   before->Add(new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
1801   after->Add( new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
1802
1803   before->Add(new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
1804   after->Add( new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
1805
1806   TH1::AddDirectory(adddirstatus);
1807 }
1808
1809 //-----------------------------------------------------------------------
1810 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
1811 {
1812   //get the number of tracks in the input event according source
1813   //selection (ESD tracks, tracklets, MC particles etc.)
1814   AliESDEvent* esd=NULL;
1815   AliAODEvent* aod=NULL;  // XZhang 20120615
1816   switch (fParamType)
1817   {
1818     case kSPDtracklet:
1819       if (!fEvent) return 0;                                           // XZhang 20120615
1820       esd = dynamic_cast<AliESDEvent*>(fEvent);
1821       aod = dynamic_cast<AliAODEvent*>(fEvent);                        // XZhang 20120615
1822 //    if (!esd) return 0;                                              // XZhang 20120615
1823 //    return esd->GetMultiplicity()->GetNumberOfTracklets();           // XZhang 20120615
1824       if (esd) return esd->GetMultiplicity()->GetNumberOfTracklets();  // XZhang 20120615
1825       if (aod) return aod->GetTracklets()->GetNumberOfTracklets();     // XZhang 20120615
1826     case kMC:
1827       if (!fMCevent) return 0;
1828       return fMCevent->GetNumberOfTracks();
1829     case kPMD:
1830       esd = dynamic_cast<AliESDEvent*>(fEvent);
1831       if (!esd) return 0;
1832       return esd->GetNumberOfPmdTracks();
1833     case kV0:
1834       return fgkNumberOfV0tracks;
1835     case kMUON:                                      // XZhang 20120604
1836       if (!fEvent) return 0;                         // XZhang 20120604
1837       esd = dynamic_cast<AliESDEvent*>(fEvent);      // XZhang 20120604
1838       if (esd) return esd->GetNumberOfMuonTracks();  // XZhang 20120604
1839       return fEvent->GetNumberOfTracks();  // if AOD // XZhang 20120604
1840     default:
1841       if (!fEvent) return 0;
1842       return fEvent->GetNumberOfTracks();
1843   }
1844   return 0;
1845 }
1846
1847 //-----------------------------------------------------------------------
1848 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
1849 {
1850   //get the input object according the data source selection:
1851   //(esd tracks, traclets, mc particles,etc...)
1852   AliESDEvent* esd=NULL;
1853   AliAODEvent* aod=NULL;  // XZhang 20120615
1854   switch (fParamType)
1855   {
1856     case kSPDtracklet:
1857       if (!fEvent) return NULL;                                              // XZhang 20120615
1858       esd = dynamic_cast<AliESDEvent*>(fEvent);
1859       aod = dynamic_cast<AliAODEvent*>(fEvent);                              // XZhang 20120615
1860 //    if (!esd) return NULL;                                                 // XZhang 20120615
1861 //    return const_cast<AliMultiplicity*>(esd->GetMultiplicity());           // XZhang 20120615
1862       if (esd) return const_cast<AliMultiplicity*>(esd->GetMultiplicity());  // XZhang 20120615
1863       if (aod) return const_cast<AliAODTracklets*>(aod->GetTracklets());     // XZhang 20120615
1864     case kMC:
1865       if (!fMCevent) return NULL;
1866       return fMCevent->GetTrack(i);
1867     case kPMD:
1868       esd = dynamic_cast<AliESDEvent*>(fEvent);
1869       if (!esd) return NULL;
1870       return esd->GetPmdTrack(i);
1871     case kV0:
1872       //esd = dynamic_cast<AliESDEvent*>(fEvent);
1873       //if (!esd) //contributed by G.Ortona
1874       //{
1875       //  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fEvent);
1876       //  if(!aod)return NULL;
1877       //  return aod->GetVZEROData();
1878       //}
1879       //return esd->GetVZEROData();
1880       return fEvent; // left only for compatibility
1881     case kMUON:                                  // XZhang 20120604
1882       if (!fEvent) return NULL;                  // XZhang 20120604
1883       esd = dynamic_cast<AliESDEvent*>(fEvent);  // XZhang 20120604
1884       if (esd) return esd->GetMuonTrack(i);      // XZhang 20120604
1885       return fEvent->GetTrack(i);  // if AOD     // XZhang 20120604
1886     default:
1887       if (!fEvent) return NULL;
1888       return fEvent->GetTrack(i);
1889   }
1890 }
1891
1892 //-----------------------------------------------------------------------
1893 void AliFlowTrackCuts::Clear(Option_t*)
1894 {
1895   //clean up
1896   fTrack=NULL;
1897   fMCevent=NULL;
1898   fMCparticle=NULL;
1899   fTrackLabel=-997;
1900   fTrackWeight=1.0;
1901   fTrackEta=0.0;
1902   fTrackPhi=0.0;
1903 }
1904
1905 //-----------------------------------------------------------------------
1906 Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
1907 {
1908   //check if passes the selected pid cut for ESDs
1909   Bool_t pass = kTRUE;
1910   switch (fPIDsource)
1911   {
1912     case kTPCpid:
1913       if (!PassesTPCpidCut(track)) pass=kFALSE;
1914       break;
1915     case kTPCdedx:
1916       if (!PassesTPCdedxCut(track)) pass=kFALSE;
1917       break;
1918     case kTOFpid:
1919       if (!PassesTOFpidCut(track)) pass=kFALSE;
1920       break;
1921     case kTOFbeta:
1922       if (!PassesTOFbetaCut(track)) pass=kFALSE;
1923       break;
1924     case kTOFbetaSimple:
1925       if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
1926       break;
1927     case kTPCbayesian:
1928       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
1929       break;
1930       // part added by F. Noferini
1931     case kTOFbayesian:
1932       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
1933       break;
1934       // end part added by F. Noferini
1935
1936       //part added by Natasha
1937     case kTPCNuclei:
1938       if (!PassesNucleiSelection(track)) pass=kFALSE;
1939       break;
1940       //end part added by Natasha
1941
1942     default:
1943       printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
1944       pass=kFALSE;
1945       break;
1946   }
1947   return pass;
1948 }
1949
1950 //-----------------------------------------------------------------------
1951 Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
1952 {
1953   //check if passes PID cut using timing in TOF
1954   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
1955                      (track->GetTOFsignal() > 12000) &&
1956                      (track->GetTOFsignal() < 100000) &&
1957                      (track->GetIntegratedLength() > 365);
1958
1959   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
1960
1961   Bool_t statusMatchingHard = TPCTOFagree(track);
1962   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
1963        return kFALSE;
1964
1965   if (!goodtrack) return kFALSE;
1966
1967   const Float_t c = 2.99792457999999984e-02;
1968   Float_t p = track->GetP();
1969   Float_t l = track->GetIntegratedLength();
1970   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
1971   Float_t timeTOF = track->GetTOFsignal()- trackT0;
1972   Float_t beta = l/timeTOF/c;
1973   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
1974   track->GetIntegratedTimes(integratedTimes);
1975   Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
1976   Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
1977   for (Int_t i=0;i<5;i++)
1978   {
1979     betaHypothesis[i] = l/integratedTimes[i]/c;
1980     s[i] = beta-betaHypothesis[i];
1981   }
1982
1983   switch (fParticleID)
1984   {
1985     case AliPID::kPion:
1986       return ( (s[2]<0.015) && (s[2]>-0.015) &&
1987                (s[3]>0.025) &&
1988                (s[4]>0.03) );
1989     case AliPID::kKaon:
1990       return ( (s[3]<0.015) && (s[3]>-0.015) &&
1991                (s[2]<-0.03) &&
1992                (s[4]>0.03) );
1993     case AliPID::kProton:
1994       return ( (s[4]<0.015) && (s[4]>-0.015) &&
1995                (s[3]<-0.025) &&
1996                (s[2]<-0.025) );
1997     default:
1998       return kFALSE;
1999   }
2000   return kFALSE;
2001 }
2002
2003 //-----------------------------------------------------------------------
2004 Float_t AliFlowTrackCuts::GetBeta(const AliESDtrack* track)
2005 {
2006   //get beta
2007   const Float_t c = 2.99792457999999984e-02;
2008   Float_t p = track->GetP();
2009   Float_t l = track->GetIntegratedLength();
2010   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2011   Float_t timeTOF = track->GetTOFsignal()- trackT0;
2012   return l/timeTOF/c;
2013 }
2014
2015 //-----------------------------------------------------------------------
2016 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
2017 {
2018   //check if track passes pid selection with an asymmetric TOF beta cut
2019   if (!fTOFpidCuts)
2020   {
2021     //printf("no TOFpidCuts\n");
2022     return kFALSE;
2023   }
2024
2025   //check if passes PID cut using timing in TOF
2026   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2027                      (track->GetTOFsignal() > 12000) &&
2028                      (track->GetTOFsignal() < 100000) &&
2029                      (track->GetIntegratedLength() > 365);
2030
2031   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2032
2033   Bool_t statusMatchingHard = TPCTOFagree(track);
2034   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2035        return kFALSE;
2036
2037   if (!goodtrack) return kFALSE;
2038
2039   Float_t beta = GetBeta(track);
2040
2041   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
2042   track->GetIntegratedTimes(integratedTimes);
2043
2044   //construct the pid index because it's not AliPID::EParticleType
2045   Int_t pid = 0;
2046   switch (fParticleID)
2047   {
2048     case AliPID::kPion:
2049       pid=2;
2050       break;
2051     case AliPID::kKaon:
2052       pid=3;
2053       break;
2054     case AliPID::kProton:
2055       pid=4;
2056       break;
2057     default:
2058       return kFALSE;
2059   }
2060
2061   //signal to cut on
2062   const Float_t c = 2.99792457999999984e-02;
2063   Float_t l = track->GetIntegratedLength();
2064   Float_t p = track->GetP();
2065   Float_t betahypothesis = l/integratedTimes[pid]/c;
2066   Float_t betadiff = beta-betahypothesis;
2067
2068   Float_t* arr = fTOFpidCuts->GetMatrixArray();
2069   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2070   if (col<0) return kFALSE;
2071   Float_t min = (*fTOFpidCuts)(1,col);
2072   Float_t max = (*fTOFpidCuts)(2,col);
2073
2074   Bool_t pass = (betadiff>min && betadiff<max);
2075
2076   return pass;
2077 }
2078
2079 //-----------------------------------------------------------------------
2080 Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
2081 {
2082   //check if passes PID cut using default TOF pid
2083   Double_t pidTOF[AliPID::kSPECIES];
2084   track->GetTOFpid(pidTOF);
2085   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
2086   return kFALSE;
2087 }
2088
2089 //-----------------------------------------------------------------------
2090 Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
2091 {
2092   //check if passes PID cut using default TPC pid
2093   Double_t pidTPC[AliPID::kSPECIES];
2094   track->GetTPCpid(pidTPC);
2095   Double_t probablity = 0.;
2096   switch (fParticleID)
2097   {
2098     case AliPID::kPion:
2099       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
2100       break;
2101     default:
2102       probablity = pidTPC[fParticleID];
2103   }
2104   if (probablity >= fParticleProbability) return kTRUE;
2105   return kFALSE;
2106 }
2107
2108 //-----------------------------------------------------------------------
2109 Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
2110 {
2111   //get TPC dedx
2112   return track->GetTPCsignal();
2113 }
2114
2115 //-----------------------------------------------------------------------
2116 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
2117 {
2118   //check if passes PID cut using dedx signal in the TPC
2119   if (!fTPCpidCuts)
2120   {
2121     //printf("no TPCpidCuts\n");
2122     return kFALSE;
2123   }
2124
2125   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
2126   if (!tpcparam) return kFALSE;
2127   Double_t p = tpcparam->GetP();
2128   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
2129   Float_t sigTPC = track->GetTPCsignal();
2130   Float_t s = (sigTPC-sigExp)/sigExp;
2131
2132   Float_t* arr = fTPCpidCuts->GetMatrixArray();
2133   Int_t arrSize = fTPCpidCuts->GetNcols();
2134   Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
2135   if (col<0) return kFALSE;
2136   Float_t min = (*fTPCpidCuts)(1,col);
2137   Float_t max = (*fTPCpidCuts)(2,col);
2138
2139   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
2140   return (s>min && s<max);
2141 }
2142
2143 //-----------------------------------------------------------------------
2144 void AliFlowTrackCuts::InitPIDcuts()
2145 {
2146   //init matrices with PID cuts
2147   TMatrixF* t = NULL;
2148   if (!fTPCpidCuts)
2149   {
2150     if (fParticleID==AliPID::kPion)
2151     {
2152       t = new TMatrixF(3,15);
2153       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.0;
2154       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.1;
2155       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.2;
2156       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.2;
2157       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
2158       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
2159       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
2160       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
2161       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
2162       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.4;  (*t)(2,9)  =  0.05;
2163       (*t)(0,10)  = 0.70;  (*t)(1,10)  = -0.4;  (*t)(2,10)  =     0;
2164       (*t)(0,11)  = 0.75;  (*t)(1,11)  = -0.4;  (*t)(2,11)  =     0;
2165       (*t)(0,12)  = 0.80;  (*t)(1,12)  = -0.4;  (*t)(2,12)  = -0.05;
2166       (*t)(0,13)  = 0.85;  (*t)(1,13)  = -0.4;  (*t)(2,13)  = -0.1;
2167       (*t)(0,14)  = 0.90;  (*t)(1,14)  = 0;     (*t)(2,14)  =     0;
2168     }
2169     else
2170     if (fParticleID==AliPID::kKaon)
2171     {
2172       t = new TMatrixF(3,12);
2173       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.2;
2174       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  = 0.2;
2175       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.2;  (*t)(2,2)  = 0.2;
2176       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.2;  (*t)(2,3)  = 0.2;
2177       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.2;
2178       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.2;
2179       (*t)(0,6)  = 0.50;  (*t)(1,6)  =-0.05;  (*t)(2,6)  = 0.2;
2180       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.1;  (*t)(2,7)  = 0.1;
2181       (*t)(0,8)  = 0.60;  (*t)(1,8)  =-0.05;  (*t)(2,8)  = 0.1;
2182       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  = 0.15;
2183       (*t)(0,10)  = 0.70;  (*t)(1,10)  = 0.05;  (*t)(2,10)  = 0.2;
2184       (*t)(0,11)  = 0.75;  (*t)(1,11)  =    0;  (*t)(2,11)  = 0;
2185     }
2186     else
2187     if (fParticleID==AliPID::kProton)
2188     {
2189       t = new TMatrixF(3,9);
2190       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.1;  (*t)(2,0)  =  0.1;
2191       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  =  0.2;
2192       (*t)(0,2)  = 0.80;  (*t)(1,2)  = -0.1;  (*t)(2,2)  =  0.2;
2193       (*t)(0,3)  = 0.85;  (*t)(1,3)  =-0.05;  (*t)(2,3)  =  0.2;
2194       (*t)(0,4)  = 0.90;  (*t)(1,4)  =-0.05;  (*t)(2,4)  = 0.25;
2195       (*t)(0,5)  = 0.95;  (*t)(1,5)  =-0.05;  (*t)(2,5)  = 0.25;
2196       (*t)(0,6)  = 1.00;  (*t)(1,6)  = -0.1;  (*t)(2,6)  = 0.25;
2197       (*t)(0,7)  = 1.10;  (*t)(1,7)  =-0.05;  (*t)(2,7)  =  0.3;
2198       (*t)(0,8) = 1.20;   (*t)(1,8)  =    0;  (*t)(2,8) =    0;
2199     }
2200     delete fTPCpidCuts;
2201     fTPCpidCuts=t;
2202   }
2203   t = NULL;
2204   if (!fTOFpidCuts)
2205   {
2206     if (fParticleID==AliPID::kPion)
2207     {
2208       //TOF pions, 0.9 purity
2209       t = new TMatrixF(3,61);
2210       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2211       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2212       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2213       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2214       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.030;  (*t)(2,4)  =   0.030;
2215       (*t)(0,5)  = 0.250;  (*t)(1,5)  = -0.036;  (*t)(2,5)  =   0.032;
2216       (*t)(0,6)  = 0.300;  (*t)(1,6)  = -0.038;  (*t)(2,6)  =   0.032;
2217       (*t)(0,7)  = 0.350;  (*t)(1,7)  = -0.034;  (*t)(2,7)  =   0.032;
2218       (*t)(0,8)  = 0.400;  (*t)(1,8)  = -0.032;  (*t)(2,8)  =   0.020;
2219       (*t)(0,9)  = 0.450;  (*t)(1,9)  = -0.030;  (*t)(2,9)  =   0.020;
2220       (*t)(0,10)  = 0.500;  (*t)(1,10)  = -0.030;  (*t)(2,10)  =   0.020;
2221       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.030;  (*t)(2,11)  =   0.020;
2222       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.030;  (*t)(2,12)  =   0.020;
2223       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.030;  (*t)(2,13)  =   0.020;
2224       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.030;  (*t)(2,14)  =   0.020;
2225       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.030;  (*t)(2,15)  =   0.020;
2226       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.030;  (*t)(2,16)  =   0.020;
2227       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.030;  (*t)(2,17)  =   0.020;
2228       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.030;  (*t)(2,18)  =   0.020;
2229       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.028;  (*t)(2,19)  =   0.028;
2230       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.028;  (*t)(2,20)  =   0.028;
2231       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.028;  (*t)(2,21)  =   0.028;
2232       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.028;
2233       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.028;
2234       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.020;  (*t)(2,24)  =   0.028;
2235       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.018;  (*t)(2,25)  =   0.028;
2236       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.016;  (*t)(2,26)  =   0.028;
2237       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.014;  (*t)(2,27)  =   0.028;
2238       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.012;  (*t)(2,28)  =   0.026;
2239       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.010;  (*t)(2,29)  =   0.026;
2240       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.008;  (*t)(2,30)  =   0.026;
2241       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.008;  (*t)(2,31)  =   0.024;
2242       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.006;  (*t)(2,32)  =   0.024;
2243       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.004;  (*t)(2,33)  =   0.024;
2244       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.004;  (*t)(2,34)  =   0.024;
2245       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.002;  (*t)(2,35)  =   0.024;
2246       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.002;  (*t)(2,36)  =   0.024;
2247       (*t)(0,37)  = 2.700;  (*t)(1,37)  = 0.000;  (*t)(2,37)  =   0.024;
2248       (*t)(0,38)  = 2.800;  (*t)(1,38)  = 0.000;  (*t)(2,38)  =   0.026;
2249       (*t)(0,39)  = 2.900;  (*t)(1,39)  = 0.000;  (*t)(2,39)  =   0.024;
2250       (*t)(0,40)  = 3.000;  (*t)(1,40)  = 0.002;  (*t)(2,40)  =   0.026;
2251       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.002;  (*t)(2,41)  =   0.026;
2252       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.002;  (*t)(2,42)  =   0.026;
2253       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.002;  (*t)(2,43)  =   0.026;
2254       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.002;  (*t)(2,44)  =   0.026;
2255       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.002;  (*t)(2,45)  =   0.026;
2256       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.002;  (*t)(2,46)  =   0.026;
2257       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.002;  (*t)(2,47)  =   0.026;
2258       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.002;  (*t)(2,48)  =   0.026;
2259       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.004;  (*t)(2,49)  =   0.024;
2260       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.004;  (*t)(2,50)  =   0.026;
2261       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.004;  (*t)(2,51)  =   0.026;
2262       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.004;  (*t)(2,52)  =   0.024;
2263       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.006;  (*t)(2,53)  =   0.024;
2264       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2265       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2266       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2267       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2268       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2269       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2270       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2271     }
2272     else
2273     if (fParticleID==AliPID::kProton)
2274     {
2275       //TOF protons, 0.9 purity
2276       t = new TMatrixF(3,61);
2277       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2278       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2279       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2280       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2281       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.07;  (*t)(2,4)  =   0.07;
2282       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.07;  (*t)(2,5)  =   0.07;
2283       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.07;  (*t)(2,6)  =   0.07;
2284       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.07;  (*t)(2,7)  =   0.07;
2285       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.07;  (*t)(2,8)  =   0.07;
2286       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.07;  (*t)(2,9)  =   0.07;
2287       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.07;  (*t)(2,10)  =   0.07;
2288       (*t)(0,11)  = 0.200;  (*t)(1,11)  = -0.07;  (*t)(2,11)  =   0.07;
2289       (*t)(0,12)  = 0.200;  (*t)(1,12)  = -0.07;  (*t)(2,12)  =   0.07;
2290       (*t)(0,13)  = 0.200;  (*t)(1,13)  = -0.07;  (*t)(2,13)  =   0.07;
2291       (*t)(0,14)  = 0.200;  (*t)(1,14)  = -0.07;  (*t)(2,14)  =   0.07;
2292       (*t)(0,15)  = 0.200;  (*t)(1,15)  = -0.07;  (*t)(2,15)  =   0.07;
2293       (*t)(0,16)  = 0.200;  (*t)(1,16)  = -0.07;  (*t)(2,16)  =   0.07;
2294       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.070;  (*t)(2,17)  =   0.070;
2295       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.072;  (*t)(2,18)  =   0.072;
2296       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.072;  (*t)(2,19)  =   0.072;
2297       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.074;  (*t)(2,20)  =   0.074;
2298       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.032;  (*t)(2,21)  =   0.032;
2299       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.026;
2300       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.026;  (*t)(2,23)  =   0.026;
2301       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.024;  (*t)(2,24)  =   0.024;
2302       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.024;  (*t)(2,25)  =   0.024;
2303       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.026;  (*t)(2,26)  =   0.026;
2304       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.026;  (*t)(2,27)  =   0.026;
2305       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.026;  (*t)(2,28)  =   0.026;
2306       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.026;  (*t)(2,29)  =   0.026;
2307       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.026;  (*t)(2,30)  =   0.026;
2308       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.026;
2309       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.026;  (*t)(2,32)  =   0.024;
2310       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.028;  (*t)(2,33)  =   0.022;
2311       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.028;  (*t)(2,34)  =   0.020;
2312       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.028;  (*t)(2,35)  =   0.018;
2313       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.028;  (*t)(2,36)  =   0.016;
2314       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.028;  (*t)(2,37)  =   0.016;
2315       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.030;  (*t)(2,38)  =   0.014;
2316       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.030;  (*t)(2,39)  =   0.012;
2317       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.030;  (*t)(2,40)  =   0.012;
2318       (*t)(0,41)  = 3.100;  (*t)(1,41)  = -0.030;  (*t)(2,41)  =   0.010;
2319       (*t)(0,42)  = 3.200;  (*t)(1,42)  = -0.030;  (*t)(2,42)  =   0.010;
2320       (*t)(0,43)  = 3.300;  (*t)(1,43)  = -0.030;  (*t)(2,43)  =   0.010;
2321       (*t)(0,44)  = 3.400;  (*t)(1,44)  = -0.030;  (*t)(2,44)  =   0.008;
2322       (*t)(0,45)  = 3.500;  (*t)(1,45)  = -0.030;  (*t)(2,45)  =   0.008;
2323       (*t)(0,46)  = 3.600;  (*t)(1,46)  = -0.030;  (*t)(2,46)  =   0.008;
2324       (*t)(0,47)  = 3.700;  (*t)(1,47)  = -0.030;  (*t)(2,47)  =   0.006;
2325       (*t)(0,48)  = 3.800;  (*t)(1,48)  = -0.030;  (*t)(2,48)  =   0.006;
2326       (*t)(0,49)  = 3.900;  (*t)(1,49)  = -0.030;  (*t)(2,49)  =   0.006;
2327       (*t)(0,50)  = 4.000;  (*t)(1,50)  = -0.028;  (*t)(2,50)  =   0.004;
2328       (*t)(0,51)  = 4.100;  (*t)(1,51)  = -0.030;  (*t)(2,51)  =   0.004;
2329       (*t)(0,52)  = 4.200;  (*t)(1,52)  = -0.030;  (*t)(2,52)  =   0.004;
2330       (*t)(0,53)  = 4.300;  (*t)(1,53)  = -0.028;  (*t)(2,53)  =   0.002;
2331       (*t)(0,54)  = 4.400;  (*t)(1,54)  = -0.030;  (*t)(2,54)  =   0.002;
2332       (*t)(0,55)  = 4.500;  (*t)(1,55)  = -0.028;  (*t)(2,55)  =   0.002;
2333       (*t)(0,56)  = 4.600;  (*t)(1,56)  = -0.028;  (*t)(2,56)  =   0.002;
2334       (*t)(0,57)  = 4.700;  (*t)(1,57)  = -0.028;  (*t)(2,57)  =   0.000;
2335       (*t)(0,58)  = 4.800;  (*t)(1,58)  = -0.028;  (*t)(2,58)  =   0.002;
2336       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2337       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2338     }
2339     else
2340     if (fParticleID==AliPID::kKaon)
2341     {
2342       //TOF kaons, 0.9 purity
2343       t = new TMatrixF(3,61);
2344       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2345       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2346       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2347       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2348       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.05;  (*t)(2,4)  =   0.05;
2349       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.05;  (*t)(2,5)  =   0.05;
2350       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.05;  (*t)(2,6)  =   0.05;
2351       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =   0.05;
2352       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.05;  (*t)(2,8)  =   0.05;
2353       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.05;  (*t)(2,9)  =   0.05;
2354       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.05;  (*t)(2,10)  =   0.05;
2355       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.026;  (*t)(2,11)  =   0.026;
2356       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.026;  (*t)(2,12)  =   0.026;
2357       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.026;  (*t)(2,13)  =   0.026;
2358       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.026;  (*t)(2,14)  =   0.026;
2359       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.026;  (*t)(2,15)  =   0.026;
2360       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.026;  (*t)(2,16)  =   0.026;
2361       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.024;  (*t)(2,17)  =   0.024;
2362       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.024;  (*t)(2,18)  =   0.024;
2363       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.024;  (*t)(2,19)  =   0.024;
2364       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.024;  (*t)(2,20)  =   0.024;
2365       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.024;  (*t)(2,21)  =   0.024;
2366       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.024;  (*t)(2,22)  =   0.022;
2367       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.020;
2368       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.026;  (*t)(2,24)  =   0.016;
2369       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.028;  (*t)(2,25)  =   0.014;
2370       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.028;  (*t)(2,26)  =   0.012;
2371       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.028;  (*t)(2,27)  =   0.010;
2372       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.028;  (*t)(2,28)  =   0.010;
2373       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.028;  (*t)(2,29)  =   0.008;
2374       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.028;  (*t)(2,30)  =   0.006;
2375       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.006;
2376       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.024;  (*t)(2,32)  =   0.004;
2377       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.020;  (*t)(2,33)  =   0.002;
2378       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.020;  (*t)(2,34)  =   0.002;
2379       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.018;  (*t)(2,35)  =   0.000;
2380       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.016;  (*t)(2,36)  =   0.000;
2381       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.014;  (*t)(2,37)  =   -0.002;
2382       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.014;  (*t)(2,38)  =   -0.004;
2383       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.012;  (*t)(2,39)  =   -0.004;
2384       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.010;  (*t)(2,40)  =   -0.006;
2385       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.000;  (*t)(2,41)  =   0.000;
2386       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.000;  (*t)(2,42)  =   0.000;
2387       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.000;  (*t)(2,43)  =   0.000;
2388       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.000;  (*t)(2,44)  =   0.000;
2389       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.000;  (*t)(2,45)  =   0.000;
2390       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.000;  (*t)(2,46)  =   0.000;
2391       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.000;  (*t)(2,47)  =   0.000;
2392       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.000;  (*t)(2,48)  =   0.000;
2393       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.000;  (*t)(2,49)  =   0.000;
2394       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.000;  (*t)(2,50)  =   0.000;
2395       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.000;  (*t)(2,51)  =   0.000;
2396       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.000;  (*t)(2,52)  =   0.000;
2397       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.000;  (*t)(2,53)  =   0.000;
2398       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2399       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2400       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2401       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2402       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2403       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2404       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2405     }
2406     delete fTOFpidCuts;
2407     fTOFpidCuts=t;
2408   }
2409 }
2410
2411 //-----------------------------------------------------------------------
2412 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
2413 {
2414   //cut on TPC bayesian pid
2415   //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
2416
2417   //Bool_t statusMatchingHard = TPCTOFagree(track);
2418   //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2419   //     return kFALSE;
2420   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
2421   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
2422   Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
2423
2424   if(! kTPC) return kFALSE;
2425
2426   Bool_t statusMatchingHard = 1;
2427   Float_t mismProb = 0;
2428   if(kTOF){
2429     statusMatchingHard = TPCTOFagree(track);
2430     mismProb = fBayesianResponse->GetTOFMismProb();
2431   }
2432   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2433        return kFALSE;
2434
2435   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
2436
2437   fProbBayes = 0.0;
2438
2439   switch (fParticleID)
2440   {
2441     case AliPID::kPion:
2442       fProbBayes = probabilities[2];
2443       break;
2444     case AliPID::kKaon:
2445       fProbBayes = probabilities[3];
2446      break;
2447     case AliPID::kProton:
2448       fProbBayes = probabilities[4];
2449       break;
2450     case AliPID::kElectron:
2451        fProbBayes = probabilities[0];
2452      break;
2453     case AliPID::kMuon:
2454        fProbBayes = probabilities[1];
2455      break;
2456     case AliPID::kDeuteron:
2457        fProbBayes = probabilities[5];
2458      break;
2459     case AliPID::kTriton:
2460        fProbBayes = probabilities[6];
2461      break;
2462     case AliPID::kHe3:
2463        fProbBayes = probabilities[7];
2464      break;
2465     default:
2466       return kFALSE;
2467   }
2468
2469   if(fProbBayes > fParticleProbability && mismProb < 0.5)
2470     {
2471       if(!fCutCharge)
2472         return kTRUE;
2473       else if (fCutCharge && fCharge * track->GetSign() > 0)
2474         return kTRUE;
2475     }
2476   return kFALSE;
2477 }
2478
2479 //-----------------------------------------------------------------------
2480 // part added by F. Noferini (some methods)
2481 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
2482 {
2483   //check is track passes bayesian combined TOF+TPC pid cut
2484   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
2485                      (track->GetStatus() & AliESDtrack::kTIME) &&
2486                      (track->GetTOFsignal() > 12000) &&
2487                      (track->GetTOFsignal() < 100000) &&
2488                      (track->GetIntegratedLength() > 365);
2489
2490   if (! goodtrack)
2491        return kFALSE;
2492
2493   if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
2494
2495   Bool_t statusMatchingHard = TPCTOFagree(track);
2496   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2497        return kFALSE;
2498
2499   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
2500   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
2501
2502   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
2503
2504   fProbBayes = 0.0;
2505
2506   switch (fParticleID)
2507   {
2508     case AliPID::kPion:
2509       fProbBayes = probabilities[2];
2510       break;
2511     case AliPID::kKaon:
2512        fProbBayes = probabilities[3];
2513      break;
2514     case AliPID::kProton:
2515        fProbBayes = probabilities[4];
2516       break;
2517     case AliPID::kElectron:
2518        fProbBayes = probabilities[0];
2519      break;
2520     case AliPID::kMuon:
2521        fProbBayes = probabilities[1];
2522      break;
2523     case AliPID::kDeuteron:
2524        fProbBayes = probabilities[5];
2525      break;
2526     case AliPID::kTriton:
2527        fProbBayes = probabilities[6];
2528      break;
2529     case AliPID::kHe3:
2530        fProbBayes = probabilities[7];
2531      break;
2532    default:
2533       return kFALSE;
2534   }
2535
2536   //  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);
2537   if(fProbBayes > fParticleProbability && mismProb < 0.5){
2538     if(!fCutCharge)
2539       return kTRUE;
2540     else if (fCutCharge && fCharge * track->GetSign() > 0)
2541       return kTRUE;
2542   }
2543   return kFALSE;
2544 }
2545
2546
2547 //-----------------------------------------------------------------------
2548  // part added by Natasha
2549 Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
2550 {
2551   //pid selection for heavy nuclei
2552   Bool_t select=kFALSE;
2553
2554   //if (!track) continue;
2555
2556   if (!track->GetInnerParam())
2557     return kFALSE;    //break;
2558
2559   const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
2560
2561   Double_t ptotTPC = tpcTrack->GetP();
2562   Double_t sigTPC = track->GetTPCsignal();
2563   Double_t dEdxBBA = 0.;
2564   Double_t dSigma = 0.;
2565
2566   switch (fParticleID)
2567   {
2568     case AliPID::kDeuteron:
2569       //pid=10;
2570       dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
2571           4.60e+00,
2572           8.9684e+00,
2573           1.640e-05,
2574           2.35e+00,
2575           2.35e+00);
2576       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
2577
2578       if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) )  )
2579       {select=kTRUE;}
2580       break;
2581
2582     case AliPID::kTriton:
2583       //pid=11;
2584       select=kFALSE;
2585       break;
2586
2587     case AliPID::kHe3:
2588       //pid=12;
2589       // ----- Pass 2 -------
2590       dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
2591           1.74962,
2592           27.4992,
2593           4.00313e-15,
2594           2.42485,
2595           8.31768);
2596       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
2597       if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
2598       {select=kTRUE;}
2599       break;
2600
2601     case AliPID::kAlpha:
2602       //pid=13;
2603       select=kFALSE;
2604       break;
2605
2606     default:
2607       return kFALSE;
2608   }
2609
2610   return select;
2611 }
2612 // end part added by Natasha
2613 //-----------------------------------------------------------------------
2614 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
2615  //set priors for the bayesian pid selection
2616   fCurrCentr = centrCur;
2617
2618   fBinLimitPID[0] = 0.300000;
2619   fBinLimitPID[1] = 0.400000;
2620   fBinLimitPID[2] = 0.500000;
2621   fBinLimitPID[3] = 0.600000;
2622   fBinLimitPID[4] = 0.700000;
2623   fBinLimitPID[5] = 0.800000;
2624   fBinLimitPID[6] = 0.900000;
2625   fBinLimitPID[7] = 1.000000;
2626   fBinLimitPID[8] = 1.200000;
2627   fBinLimitPID[9] = 1.400000;
2628   fBinLimitPID[10] = 1.600000;
2629   fBinLimitPID[11] = 1.800000;
2630   fBinLimitPID[12] = 2.000000;
2631   fBinLimitPID[13] = 2.200000;
2632   fBinLimitPID[14] = 2.400000;
2633   fBinLimitPID[15] = 2.600000;
2634   fBinLimitPID[16] = 2.800000;
2635   fBinLimitPID[17] = 3.000000;
2636
2637   // 0-10%
2638   if(centrCur < 10){
2639       fC[0][0] = 0.005;
2640       fC[0][1] = 0.005;
2641       fC[0][2] = 1.0000;
2642       fC[0][3] = 0.0010;
2643       fC[0][4] = 0.0010;
2644
2645       fC[1][0] = 0.005;
2646       fC[1][1] = 0.005;
2647       fC[1][2] = 1.0000;
2648       fC[1][3] = 0.0168;
2649       fC[1][4] = 0.0122;
2650
2651       fC[2][0] = 0.005;
2652       fC[2][1] = 0.005;
2653       fC[2][2] = 1.0000;
2654       fC[2][3] = 0.0272;
2655       fC[2][4] = 0.0070;
2656
2657       fC[3][0] = 0.005;
2658       fC[3][1] = 0.005;
2659       fC[3][2] = 1.0000;
2660       fC[3][3] = 0.0562;
2661       fC[3][4] = 0.0258;
2662
2663       fC[4][0] = 0.005;
2664       fC[4][1] = 0.005;
2665       fC[4][2] = 1.0000;
2666       fC[4][3] = 0.0861;
2667       fC[4][4] = 0.0496;
2668
2669       fC[5][0] = 0.005;
2670       fC[5][1] = 0.005;
2671       fC[5][2] = 1.0000;
2672       fC[5][3] = 0.1168;
2673       fC[5][4] = 0.0740;
2674
2675       fC[6][0] = 0.005;
2676       fC[6][1] = 0.005;
2677       fC[6][2] = 1.0000;
2678       fC[6][3] = 0.1476;
2679       fC[6][4] = 0.0998;
2680
2681       fC[7][0] = 0.005;
2682       fC[7][1] = 0.005;
2683       fC[7][2] = 1.0000;
2684       fC[7][3] = 0.1810;
2685       fC[7][4] = 0.1296;
2686
2687       fC[8][0] = 0.005;
2688       fC[8][1] = 0.005;
2689       fC[8][2] = 1.0000;
2690       fC[8][3] = 0.2240;
2691       fC[8][4] = 0.1827;
2692
2693       fC[9][0] = 0.005;
2694       fC[9][1] = 0.005;
2695       fC[9][2] = 1.0000;
2696       fC[9][3] = 0.2812;
2697       fC[9][4] = 0.2699;
2698
2699       fC[10][0] = 0.005;
2700       fC[10][1] = 0.005;
2701       fC[10][2] = 1.0000;
2702       fC[10][3] = 0.3328;
2703       fC[10][4] = 0.3714;
2704
2705       fC[11][0] = 0.005;
2706       fC[11][1] = 0.005;
2707       fC[11][2] = 1.0000;
2708       fC[11][3] = 0.3780;
2709       fC[11][4] = 0.4810;
2710
2711       fC[12][0] = 0.005;
2712       fC[12][1] = 0.005;
2713       fC[12][2] = 1.0000;
2714       fC[12][3] = 0.4125;
2715       fC[12][4] = 0.5771;
2716
2717       fC[13][0] = 0.005;
2718       fC[13][1] = 0.005;
2719       fC[13][2] = 1.0000;
2720       fC[13][3] = 0.4486;
2721       fC[13][4] = 0.6799;
2722
2723       fC[14][0] = 0.005;
2724       fC[14][1] = 0.005;
2725       fC[14][2] = 1.0000;
2726       fC[14][3] = 0.4840;
2727       fC[14][4] = 0.7668;
2728
2729       fC[15][0] = 0.005;
2730       fC[15][1] = 0.005;
2731       fC[15][2] = 1.0000;
2732       fC[15][3] = 0.4971;
2733       fC[15][4] = 0.8288;
2734
2735       fC[16][0] = 0.005;
2736       fC[16][1] = 0.005;
2737       fC[16][2] = 1.0000;
2738       fC[16][3] = 0.4956;
2739       fC[16][4] = 0.8653;
2740
2741       fC[17][0] = 0.005;
2742       fC[17][1] = 0.005;
2743       fC[17][2] = 1.0000;
2744       fC[17][3] = 0.5173;
2745       fC[17][4] = 0.9059;
2746   }
2747   // 10-20%
2748   else if(centrCur < 20){
2749      fC[0][0] = 0.005;
2750       fC[0][1] = 0.005;
2751       fC[0][2] = 1.0000;
2752       fC[0][3] = 0.0010;
2753       fC[0][4] = 0.0010;
2754
2755       fC[1][0] = 0.005;
2756       fC[1][1] = 0.005;
2757       fC[1][2] = 1.0000;
2758       fC[1][3] = 0.0132;
2759       fC[1][4] = 0.0088;
2760
2761       fC[2][0] = 0.005;
2762       fC[2][1] = 0.005;
2763       fC[2][2] = 1.0000;
2764       fC[2][3] = 0.0283;
2765       fC[2][4] = 0.0068;
2766
2767       fC[3][0] = 0.005;
2768       fC[3][1] = 0.005;
2769       fC[3][2] = 1.0000;
2770       fC[3][3] = 0.0577;
2771       fC[3][4] = 0.0279;
2772
2773       fC[4][0] = 0.005;
2774       fC[4][1] = 0.005;
2775       fC[4][2] = 1.0000;
2776       fC[4][3] = 0.0884;
2777       fC[4][4] = 0.0534;
2778
2779       fC[5][0] = 0.005;
2780       fC[5][1] = 0.005;
2781       fC[5][2] = 1.0000;
2782       fC[5][3] = 0.1179;
2783       fC[5][4] = 0.0794;
2784
2785       fC[6][0] = 0.005;
2786       fC[6][1] = 0.005;
2787       fC[6][2] = 1.0000;
2788       fC[6][3] = 0.1480;
2789       fC[6][4] = 0.1058;
2790
2791       fC[7][0] = 0.005;
2792       fC[7][1] = 0.005;
2793       fC[7][2] = 1.0000;
2794       fC[7][3] = 0.1807;
2795       fC[7][4] = 0.1366;
2796
2797       fC[8][0] = 0.005;
2798       fC[8][1] = 0.005;
2799       fC[8][2] = 1.0000;
2800       fC[8][3] = 0.2219;
2801       fC[8][4] = 0.1891;
2802
2803       fC[9][0] = 0.005;
2804       fC[9][1] = 0.005;
2805       fC[9][2] = 1.0000;
2806       fC[9][3] = 0.2804;
2807       fC[9][4] = 0.2730;
2808
2809       fC[10][0] = 0.005;
2810       fC[10][1] = 0.005;
2811       fC[10][2] = 1.0000;
2812       fC[10][3] = 0.3283;
2813       fC[10][4] = 0.3660;
2814
2815       fC[11][0] = 0.005;
2816       fC[11][1] = 0.005;
2817       fC[11][2] = 1.0000;
2818       fC[11][3] = 0.3710;
2819       fC[11][4] = 0.4647;
2820
2821       fC[12][0] = 0.005;
2822       fC[12][1] = 0.005;
2823       fC[12][2] = 1.0000;
2824       fC[12][3] = 0.4093;
2825       fC[12][4] = 0.5566;
2826
2827       fC[13][0] = 0.005;
2828       fC[13][1] = 0.005;
2829       fC[13][2] = 1.0000;
2830       fC[13][3] = 0.4302;
2831       fC[13][4] = 0.6410;
2832
2833       fC[14][0] = 0.005;
2834       fC[14][1] = 0.005;
2835       fC[14][2] = 1.0000;
2836       fC[14][3] = 0.4649;
2837       fC[14][4] = 0.7055;
2838
2839       fC[15][0] = 0.005;
2840       fC[15][1] = 0.005;
2841       fC[15][2] = 1.0000;
2842       fC[15][3] = 0.4523;
2843       fC[15][4] = 0.7440;
2844
2845       fC[16][0] = 0.005;
2846       fC[16][1] = 0.005;
2847       fC[16][2] = 1.0000;
2848       fC[16][3] = 0.4591;
2849       fC[16][4] = 0.7799;
2850
2851       fC[17][0] = 0.005;
2852       fC[17][1] = 0.005;
2853       fC[17][2] = 1.0000;
2854       fC[17][3] = 0.4804;
2855       fC[17][4] = 0.8218;
2856   }
2857   // 20-30%
2858   else if(centrCur < 30){
2859      fC[0][0] = 0.005;
2860       fC[0][1] = 0.005;
2861       fC[0][2] = 1.0000;
2862       fC[0][3] = 0.0010;
2863       fC[0][4] = 0.0010;
2864
2865       fC[1][0] = 0.005;
2866       fC[1][1] = 0.005;
2867       fC[1][2] = 1.0000;
2868       fC[1][3] = 0.0102;
2869       fC[1][4] = 0.0064;
2870
2871       fC[2][0] = 0.005;
2872       fC[2][1] = 0.005;
2873       fC[2][2] = 1.0000;
2874       fC[2][3] = 0.0292;
2875       fC[2][4] = 0.0066;
2876
2877       fC[3][0] = 0.005;
2878       fC[3][1] = 0.005;
2879       fC[3][2] = 1.0000;
2880       fC[3][3] = 0.0597;
2881       fC[3][4] = 0.0296;
2882
2883       fC[4][0] = 0.005;
2884       fC[4][1] = 0.005;
2885       fC[4][2] = 1.0000;
2886       fC[4][3] = 0.0900;
2887       fC[4][4] = 0.0589;
2888
2889       fC[5][0] = 0.005;
2890       fC[5][1] = 0.005;
2891       fC[5][2] = 1.0000;
2892       fC[5][3] = 0.1199;
2893       fC[5][4] = 0.0859;
2894
2895       fC[6][0] = 0.005;
2896       fC[6][1] = 0.005;
2897       fC[6][2] = 1.0000;
2898       fC[6][3] = 0.1505;
2899       fC[6][4] = 0.1141;
2900
2901       fC[7][0] = 0.005;
2902       fC[7][1] = 0.005;
2903       fC[7][2] = 1.0000;
2904       fC[7][3] = 0.1805;
2905       fC[7][4] = 0.1454;
2906
2907       fC[8][0] = 0.005;
2908       fC[8][1] = 0.005;
2909       fC[8][2] = 1.0000;
2910       fC[8][3] = 0.2221;
2911       fC[8][4] = 0.2004;
2912
2913       fC[9][0] = 0.005;
2914       fC[9][1] = 0.005;
2915       fC[9][2] = 1.0000;
2916       fC[9][3] = 0.2796;
2917       fC[9][4] = 0.2838;
2918
2919       fC[10][0] = 0.005;
2920       fC[10][1] = 0.005;
2921       fC[10][2] = 1.0000;
2922       fC[10][3] = 0.3271;
2923       fC[10][4] = 0.3682;
2924
2925       fC[11][0] = 0.005;
2926       fC[11][1] = 0.005;
2927       fC[11][2] = 1.0000;
2928       fC[11][3] = 0.3648;
2929       fC[11][4] = 0.4509;
2930
2931       fC[12][0] = 0.005;
2932       fC[12][1] = 0.005;
2933       fC[12][2] = 1.0000;
2934       fC[12][3] = 0.3988;
2935       fC[12][4] = 0.5339;
2936
2937       fC[13][0] = 0.005;
2938       fC[13][1] = 0.005;
2939       fC[13][2] = 1.0000;
2940       fC[13][3] = 0.4315;
2941       fC[13][4] = 0.5995;
2942
2943       fC[14][0] = 0.005;
2944       fC[14][1] = 0.005;
2945       fC[14][2] = 1.0000;
2946       fC[14][3] = 0.4548;
2947       fC[14][4] = 0.6612;
2948
2949       fC[15][0] = 0.005;
2950       fC[15][1] = 0.005;
2951       fC[15][2] = 1.0000;
2952       fC[15][3] = 0.4744;
2953       fC[15][4] = 0.7060;
2954
2955       fC[16][0] = 0.005;
2956       fC[16][1] = 0.005;
2957       fC[16][2] = 1.0000;
2958       fC[16][3] = 0.4899;
2959       fC[16][4] = 0.7388;
2960
2961       fC[17][0] = 0.005;
2962       fC[17][1] = 0.005;
2963       fC[17][2] = 1.0000;
2964       fC[17][3] = 0.4411;
2965       fC[17][4] = 0.7293;
2966   }
2967   // 30-40%
2968   else if(centrCur < 40){
2969       fC[0][0] = 0.005;
2970       fC[0][1] = 0.005;
2971       fC[0][2] = 1.0000;
2972       fC[0][3] = 0.0010;
2973       fC[0][4] = 0.0010;
2974
2975       fC[1][0] = 0.005;
2976       fC[1][1] = 0.005;
2977       fC[1][2] = 1.0000;
2978       fC[1][3] = 0.0102;
2979       fC[1][4] = 0.0048;
2980
2981       fC[2][0] = 0.005;
2982       fC[2][1] = 0.005;
2983       fC[2][2] = 1.0000;
2984       fC[2][3] = 0.0306;
2985       fC[2][4] = 0.0079;
2986
2987       fC[3][0] = 0.005;
2988       fC[3][1] = 0.005;
2989       fC[3][2] = 1.0000;
2990       fC[3][3] = 0.0617;
2991       fC[3][4] = 0.0338;
2992
2993       fC[4][0] = 0.005;
2994       fC[4][1] = 0.005;
2995       fC[4][2] = 1.0000;
2996       fC[4][3] = 0.0920;
2997       fC[4][4] = 0.0652;
2998
2999       fC[5][0] = 0.005;
3000       fC[5][1] = 0.005;
3001       fC[5][2] = 1.0000;
3002       fC[5][3] = 0.1211;
3003       fC[5][4] = 0.0955;
3004
3005       fC[6][0] = 0.005;
3006       fC[6][1] = 0.005;
3007       fC[6][2] = 1.0000;
3008       fC[6][3] = 0.1496;
3009       fC[6][4] = 0.1242;
3010
3011       fC[7][0] = 0.005;
3012       fC[7][1] = 0.005;
3013       fC[7][2] = 1.0000;
3014       fC[7][3] = 0.1807;
3015       fC[7][4] = 0.1576;
3016
3017       fC[8][0] = 0.005;
3018       fC[8][1] = 0.005;
3019       fC[8][2] = 1.0000;
3020       fC[8][3] = 0.2195;
3021       fC[8][4] = 0.2097;
3022
3023       fC[9][0] = 0.005;
3024       fC[9][1] = 0.005;
3025       fC[9][2] = 1.0000;
3026       fC[9][3] = 0.2732;
3027       fC[9][4] = 0.2884;
3028
3029       fC[10][0] = 0.005;
3030       fC[10][1] = 0.005;
3031       fC[10][2] = 1.0000;
3032       fC[10][3] = 0.3204;
3033       fC[10][4] = 0.3679;
3034
3035       fC[11][0] = 0.005;
3036       fC[11][1] = 0.005;
3037       fC[11][2] = 1.0000;
3038       fC[11][3] = 0.3564;
3039       fC[11][4] = 0.4449;
3040
3041       fC[12][0] = 0.005;
3042       fC[12][1] = 0.005;
3043       fC[12][2] = 1.0000;
3044       fC[12][3] = 0.3791;
3045       fC[12][4] = 0.5052;
3046
3047       fC[13][0] = 0.005;
3048       fC[13][1] = 0.005;
3049       fC[13][2] = 1.0000;
3050       fC[13][3] = 0.4062;
3051       fC[13][4] = 0.5647;
3052
3053       fC[14][0] = 0.005;
3054       fC[14][1] = 0.005;
3055       fC[14][2] = 1.0000;
3056       fC[14][3] = 0.4234;
3057       fC[14][4] = 0.6203;
3058
3059       fC[15][0] = 0.005;
3060       fC[15][1] = 0.005;
3061       fC[15][2] = 1.0000;
3062       fC[15][3] = 0.4441;
3063       fC[15][4] = 0.6381;
3064
3065       fC[16][0] = 0.005;
3066       fC[16][1] = 0.005;
3067       fC[16][2] = 1.0000;
3068       fC[16][3] = 0.4629;
3069       fC[16][4] = 0.6496;
3070
3071       fC[17][0] = 0.005;
3072       fC[17][1] = 0.005;
3073       fC[17][2] = 1.0000;
3074       fC[17][3] = 0.4293;
3075       fC[17][4] = 0.6491;
3076   }
3077   // 40-50%
3078   else if(centrCur < 50){
3079       fC[0][0] = 0.005;
3080       fC[0][1] = 0.005;
3081       fC[0][2] = 1.0000;
3082       fC[0][3] = 0.0010;
3083       fC[0][4] = 0.0010;
3084
3085       fC[1][0] = 0.005;
3086       fC[1][1] = 0.005;
3087       fC[1][2] = 1.0000;
3088       fC[1][3] = 0.0093;
3089       fC[1][4] = 0.0057;
3090
3091       fC[2][0] = 0.005;
3092       fC[2][1] = 0.005;
3093       fC[2][2] = 1.0000;
3094       fC[2][3] = 0.0319;
3095       fC[2][4] = 0.0075;
3096
3097       fC[3][0] = 0.005;
3098       fC[3][1] = 0.005;
3099       fC[3][2] = 1.0000;
3100       fC[3][3] = 0.0639;
3101       fC[3][4] = 0.0371;
3102
3103       fC[4][0] = 0.005;
3104       fC[4][1] = 0.005;
3105       fC[4][2] = 1.0000;
3106       fC[4][3] = 0.0939;
3107       fC[4][4] = 0.0725;
3108
3109       fC[5][0] = 0.005;
3110       fC[5][1] = 0.005;
3111       fC[5][2] = 1.0000;
3112       fC[5][3] = 0.1224;
3113       fC[5][4] = 0.1045;
3114
3115       fC[6][0] = 0.005;
3116       fC[6][1] = 0.005;
3117       fC[6][2] = 1.0000;
3118       fC[6][3] = 0.1520;
3119       fC[6][4] = 0.1387;
3120
3121       fC[7][0] = 0.005;
3122       fC[7][1] = 0.005;
3123       fC[7][2] = 1.0000;
3124       fC[7][3] = 0.1783;
3125       fC[7][4] = 0.1711;
3126
3127       fC[8][0] = 0.005;
3128       fC[8][1] = 0.005;
3129       fC[8][2] = 1.0000;
3130       fC[8][3] = 0.2202;
3131       fC[8][4] = 0.2269;
3132
3133       fC[9][0] = 0.005;
3134       fC[9][1] = 0.005;
3135       fC[9][2] = 1.0000;
3136       fC[9][3] = 0.2672;
3137       fC[9][4] = 0.2955;
3138
3139       fC[10][0] = 0.005;
3140       fC[10][1] = 0.005;
3141       fC[10][2] = 1.0000;
3142       fC[10][3] = 0.3191;
3143       fC[10][4] = 0.3676;
3144
3145       fC[11][0] = 0.005;
3146       fC[11][1] = 0.005;
3147       fC[11][2] = 1.0000;
3148       fC[11][3] = 0.3434;
3149       fC[11][4] = 0.4321;
3150
3151       fC[12][0] = 0.005;
3152       fC[12][1] = 0.005;
3153       fC[12][2] = 1.0000;
3154       fC[12][3] = 0.3692;
3155       fC[12][4] = 0.4879;
3156
3157       fC[13][0] = 0.005;
3158       fC[13][1] = 0.005;
3159       fC[13][2] = 1.0000;
3160       fC[13][3] = 0.3993;
3161       fC[13][4] = 0.5377;
3162
3163       fC[14][0] = 0.005;
3164       fC[14][1] = 0.005;
3165       fC[14][2] = 1.0000;
3166       fC[14][3] = 0.3818;
3167       fC[14][4] = 0.5547;
3168
3169       fC[15][0] = 0.005;
3170       fC[15][1] = 0.005;
3171       fC[15][2] = 1.0000;
3172       fC[15][3] = 0.4003;
3173       fC[15][4] = 0.5484;
3174
3175       fC[16][0] = 0.005;
3176       fC[16][1] = 0.005;
3177       fC[16][2] = 1.0000;
3178       fC[16][3] = 0.4281;
3179       fC[16][4] = 0.5383;
3180
3181       fC[17][0] = 0.005;
3182       fC[17][1] = 0.005;
3183       fC[17][2] = 1.0000;
3184       fC[17][3] = 0.3960;
3185       fC[17][4] = 0.5374;
3186   }
3187   // 50-60%
3188   else if(centrCur < 60){
3189       fC[0][0] = 0.005;
3190       fC[0][1] = 0.005;
3191       fC[0][2] = 1.0000;
3192       fC[0][3] = 0.0010;
3193       fC[0][4] = 0.0010;
3194
3195       fC[1][0] = 0.005;
3196       fC[1][1] = 0.005;
3197       fC[1][2] = 1.0000;
3198       fC[1][3] = 0.0076;
3199       fC[1][4] = 0.0032;
3200
3201       fC[2][0] = 0.005;
3202       fC[2][1] = 0.005;
3203       fC[2][2] = 1.0000;
3204       fC[2][3] = 0.0329;
3205       fC[2][4] = 0.0085;
3206
3207       fC[3][0] = 0.005;
3208       fC[3][1] = 0.005;
3209       fC[3][2] = 1.0000;
3210       fC[3][3] = 0.0653;
3211       fC[3][4] = 0.0423;
3212
3213       fC[4][0] = 0.005;
3214       fC[4][1] = 0.005;
3215       fC[4][2] = 1.0000;
3216       fC[4][3] = 0.0923;
3217       fC[4][4] = 0.0813;
3218
3219       fC[5][0] = 0.005;
3220       fC[5][1] = 0.005;
3221       fC[5][2] = 1.0000;
3222       fC[5][3] = 0.1219;
3223       fC[5][4] = 0.1161;
3224
3225       fC[6][0] = 0.005;
3226       fC[6][1] = 0.005;
3227       fC[6][2] = 1.0000;
3228       fC[6][3] = 0.1519;
3229       fC[6][4] = 0.1520;
3230
3231       fC[7][0] = 0.005;
3232       fC[7][1] = 0.005;
3233       fC[7][2] = 1.0000;
3234       fC[7][3] = 0.1763;
3235       fC[7][4] = 0.1858;
3236
3237       fC[8][0] = 0.005;
3238       fC[8][1] = 0.005;
3239       fC[8][2] = 1.0000;
3240       fC[8][3] = 0.2178;
3241       fC[8][4] = 0.2385;
3242
3243       fC[9][0] = 0.005;
3244       fC[9][1] = 0.005;
3245       fC[9][2] = 1.0000;
3246       fC[9][3] = 0.2618;
3247       fC[9][4] = 0.3070;
3248
3249       fC[10][0] = 0.005;
3250       fC[10][1] = 0.005;
3251       fC[10][2] = 1.0000;
3252       fC[10][3] = 0.3067;
3253       fC[10][4] = 0.3625;
3254
3255       fC[11][0] = 0.005;
3256       fC[11][1] = 0.005;
3257       fC[11][2] = 1.0000;
3258       fC[11][3] = 0.3336;
3259       fC[11][4] = 0.4188;
3260
3261       fC[12][0] = 0.005;
3262       fC[12][1] = 0.005;
3263       fC[12][2] = 1.0000;
3264       fC[12][3] = 0.3706;
3265       fC[12][4] = 0.4511;
3266
3267       fC[13][0] = 0.005;
3268       fC[13][1] = 0.005;
3269       fC[13][2] = 1.0000;
3270       fC[13][3] = 0.3765;
3271       fC[13][4] = 0.4729;
3272
3273       fC[14][0] = 0.005;
3274       fC[14][1] = 0.005;
3275       fC[14][2] = 1.0000;
3276       fC[14][3] = 0.3942;
3277       fC[14][4] = 0.4855;
3278
3279       fC[15][0] = 0.005;
3280       fC[15][1] = 0.005;
3281       fC[15][2] = 1.0000;
3282       fC[15][3] = 0.4051;
3283       fC[15][4] = 0.4762;
3284
3285       fC[16][0] = 0.005;
3286       fC[16][1] = 0.005;
3287       fC[16][2] = 1.0000;
3288       fC[16][3] = 0.3843;
3289       fC[16][4] = 0.4763;
3290
3291       fC[17][0] = 0.005;
3292       fC[17][1] = 0.005;
3293       fC[17][2] = 1.0000;
3294       fC[17][3] = 0.4237;
3295       fC[17][4] = 0.4773;
3296   }
3297   // 60-70%
3298   else if(centrCur < 70){
3299          fC[0][0] = 0.005;
3300       fC[0][1] = 0.005;
3301       fC[0][2] = 1.0000;
3302       fC[0][3] = 0.0010;
3303       fC[0][4] = 0.0010;
3304
3305       fC[1][0] = 0.005;
3306       fC[1][1] = 0.005;
3307       fC[1][2] = 1.0000;
3308       fC[1][3] = 0.0071;
3309       fC[1][4] = 0.0012;
3310
3311       fC[2][0] = 0.005;
3312       fC[2][1] = 0.005;
3313       fC[2][2] = 1.0000;
3314       fC[2][3] = 0.0336;
3315       fC[2][4] = 0.0097;
3316
3317       fC[3][0] = 0.005;
3318       fC[3][1] = 0.005;
3319       fC[3][2] = 1.0000;
3320       fC[3][3] = 0.0662;
3321       fC[3][4] = 0.0460;
3322
3323       fC[4][0] = 0.005;
3324       fC[4][1] = 0.005;
3325       fC[4][2] = 1.0000;
3326       fC[4][3] = 0.0954;
3327       fC[4][4] = 0.0902;
3328
3329       fC[5][0] = 0.005;
3330       fC[5][1] = 0.005;
3331       fC[5][2] = 1.0000;
3332       fC[5][3] = 0.1181;
3333       fC[5][4] = 0.1306;
3334
3335       fC[6][0] = 0.005;
3336       fC[6][1] = 0.005;
3337       fC[6][2] = 1.0000;
3338       fC[6][3] = 0.1481;
3339       fC[6][4] = 0.1662;
3340
3341       fC[7][0] = 0.005;
3342       fC[7][1] = 0.005;
3343       fC[7][2] = 1.0000;
3344       fC[7][3] = 0.1765;
3345       fC[7][4] = 0.1963;
3346
3347       fC[8][0] = 0.005;
3348       fC[8][1] = 0.005;
3349       fC[8][2] = 1.0000;
3350       fC[8][3] = 0.2155;
3351       fC[8][4] = 0.2433;
3352
3353       fC[9][0] = 0.005;
3354       fC[9][1] = 0.005;
3355       fC[9][2] = 1.0000;
3356       fC[9][3] = 0.2580;
3357       fC[9][4] = 0.3022;
3358
3359       fC[10][0] = 0.005;
3360       fC[10][1] = 0.005;
3361       fC[10][2] = 1.0000;
3362       fC[10][3] = 0.2872;
3363       fC[10][4] = 0.3481;
3364
3365       fC[11][0] = 0.005;
3366       fC[11][1] = 0.005;
3367       fC[11][2] = 1.0000;
3368       fC[11][3] = 0.3170;
3369       fC[11][4] = 0.3847;
3370
3371       fC[12][0] = 0.005;
3372       fC[12][1] = 0.005;
3373       fC[12][2] = 1.0000;
3374       fC[12][3] = 0.3454;
3375       fC[12][4] = 0.4258;
3376
3377       fC[13][0] = 0.005;
3378       fC[13][1] = 0.005;
3379       fC[13][2] = 1.0000;
3380       fC[13][3] = 0.3580;
3381       fC[13][4] = 0.4299;
3382
3383       fC[14][0] = 0.005;
3384       fC[14][1] = 0.005;
3385       fC[14][2] = 1.0000;
3386       fC[14][3] = 0.3903;
3387       fC[14][4] = 0.4326;
3388
3389       fC[15][0] = 0.005;
3390       fC[15][1] = 0.005;
3391       fC[15][2] = 1.0000;
3392       fC[15][3] = 0.3690;
3393       fC[15][4] = 0.4491;
3394
3395       fC[16][0] = 0.005;
3396       fC[16][1] = 0.005;
3397       fC[16][2] = 1.0000;
3398       fC[16][3] = 0.4716;
3399       fC[16][4] = 0.4298;
3400
3401       fC[17][0] = 0.005;
3402       fC[17][1] = 0.005;
3403       fC[17][2] = 1.0000;
3404       fC[17][3] = 0.3875;
3405       fC[17][4] = 0.4083;
3406   }
3407   // 70-80%
3408   else if(centrCur < 80){
3409       fC[0][0] = 0.005;
3410       fC[0][1] = 0.005;
3411       fC[0][2] = 1.0000;
3412       fC[0][3] = 0.0010;
3413       fC[0][4] = 0.0010;
3414
3415       fC[1][0] = 0.005;
3416       fC[1][1] = 0.005;
3417       fC[1][2] = 1.0000;
3418       fC[1][3] = 0.0075;
3419       fC[1][4] = 0.0007;
3420
3421       fC[2][0] = 0.005;
3422       fC[2][1] = 0.005;
3423       fC[2][2] = 1.0000;
3424       fC[2][3] = 0.0313;
3425       fC[2][4] = 0.0124;
3426
3427       fC[3][0] = 0.005;
3428       fC[3][1] = 0.005;
3429       fC[3][2] = 1.0000;
3430       fC[3][3] = 0.0640;
3431       fC[3][4] = 0.0539;
3432
3433       fC[4][0] = 0.005;
3434       fC[4][1] = 0.005;
3435       fC[4][2] = 1.0000;
3436       fC[4][3] = 0.0923;
3437       fC[4][4] = 0.0992;
3438
3439       fC[5][0] = 0.005;
3440       fC[5][1] = 0.005;
3441       fC[5][2] = 1.0000;
3442       fC[5][3] = 0.1202;
3443       fC[5][4] = 0.1417;
3444
3445       fC[6][0] = 0.005;
3446       fC[6][1] = 0.005;
3447       fC[6][2] = 1.0000;
3448       fC[6][3] = 0.1413;
3449       fC[6][4] = 0.1729;
3450
3451       fC[7][0] = 0.005;
3452       fC[7][1] = 0.005;
3453       fC[7][2] = 1.0000;
3454       fC[7][3] = 0.1705;
3455       fC[7][4] = 0.1999;
3456
3457       fC[8][0] = 0.005;
3458       fC[8][1] = 0.005;
3459       fC[8][2] = 1.0000;
3460       fC[8][3] = 0.2103;
3461       fC[8][4] = 0.2472;
3462
3463       fC[9][0] = 0.005;
3464       fC[9][1] = 0.005;
3465       fC[9][2] = 1.0000;
3466       fC[9][3] = 0.2373;
3467       fC[9][4] = 0.2916;
3468
3469       fC[10][0] = 0.005;
3470       fC[10][1] = 0.005;
3471       fC[10][2] = 1.0000;
3472       fC[10][3] = 0.2824;
3473       fC[10][4] = 0.3323;
3474
3475       fC[11][0] = 0.005;
3476       fC[11][1] = 0.005;
3477       fC[11][2] = 1.0000;
3478       fC[11][3] = 0.3046;
3479       fC[11][4] = 0.3576;
3480
3481       fC[12][0] = 0.005;
3482       fC[12][1] = 0.005;
3483       fC[12][2] = 1.0000;
3484       fC[12][3] = 0.3585;
3485       fC[12][4] = 0.4003;
3486
3487       fC[13][0] = 0.005;
3488       fC[13][1] = 0.005;
3489       fC[13][2] = 1.0000;
3490       fC[13][3] = 0.3461;
3491       fC[13][4] = 0.3982;
3492
3493       fC[14][0] = 0.005;
3494       fC[14][1] = 0.005;
3495       fC[14][2] = 1.0000;
3496       fC[14][3] = 0.3362;
3497       fC[14][4] = 0.3776;
3498
3499       fC[15][0] = 0.005;
3500       fC[15][1] = 0.005;
3501       fC[15][2] = 1.0000;
3502       fC[15][3] = 0.3071;
3503       fC[15][4] = 0.3500;
3504
3505       fC[16][0] = 0.005;
3506       fC[16][1] = 0.005;
3507       fC[16][2] = 1.0000;
3508       fC[16][3] = 0.2914;
3509       fC[16][4] = 0.3937;
3510
3511       fC[17][0] = 0.005;
3512       fC[17][1] = 0.005;
3513       fC[17][2] = 1.0000;
3514       fC[17][3] = 0.3727;
3515       fC[17][4] = 0.3877;
3516   }
3517   // 80-100%
3518   else{
3519       fC[0][0] = 0.005;
3520       fC[0][1] = 0.005;
3521       fC[0][2] = 1.0000;
3522       fC[0][3] = 0.0010;
3523       fC[0][4] = 0.0010;
3524
3525       fC[1][0] = 0.005;
3526       fC[1][1] = 0.005;
3527       fC[1][2] = 1.0000;
3528       fC[1][3] = 0.0060;
3529       fC[1][4] = 0.0035;
3530
3531       fC[2][0] = 0.005;
3532       fC[2][1] = 0.005;
3533       fC[2][2] = 1.0000;
3534       fC[2][3] = 0.0323;
3535       fC[2][4] = 0.0113;
3536
3537       fC[3][0] = 0.005;
3538       fC[3][1] = 0.005;
3539       fC[3][2] = 1.0000;
3540       fC[3][3] = 0.0609;
3541       fC[3][4] = 0.0653;
3542
3543       fC[4][0] = 0.005;
3544       fC[4][1] = 0.005;
3545       fC[4][2] = 1.0000;
3546       fC[4][3] = 0.0922;
3547       fC[4][4] = 0.1076;
3548
3549       fC[5][0] = 0.005;
3550       fC[5][1] = 0.005;
3551       fC[5][2] = 1.0000;
3552       fC[5][3] = 0.1096;
3553       fC[5][4] = 0.1328;
3554
3555       fC[6][0] = 0.005;
3556       fC[6][1] = 0.005;
3557       fC[6][2] = 1.0000;
3558       fC[6][3] = 0.1495;
3559       fC[6][4] = 0.1779;
3560
3561       fC[7][0] = 0.005;
3562       fC[7][1] = 0.005;
3563       fC[7][2] = 1.0000;
3564       fC[7][3] = 0.1519;
3565       fC[7][4] = 0.1989;
3566
3567       fC[8][0] = 0.005;
3568       fC[8][1] = 0.005;
3569       fC[8][2] = 1.0000;
3570       fC[8][3] = 0.1817;
3571       fC[8][4] = 0.2472;
3572
3573       fC[9][0] = 0.005;
3574       fC[9][1] = 0.005;
3575       fC[9][2] = 1.0000;
3576       fC[9][3] = 0.2429;
3577       fC[9][4] = 0.2684;
3578
3579       fC[10][0] = 0.005;
3580       fC[10][1] = 0.005;
3581       fC[10][2] = 1.0000;
3582       fC[10][3] = 0.2760;
3583       fC[10][4] = 0.3098;
3584
3585       fC[11][0] = 0.005;
3586       fC[11][1] = 0.005;
3587       fC[11][2] = 1.0000;
3588       fC[11][3] = 0.2673;
3589       fC[11][4] = 0.3198;
3590
3591       fC[12][0] = 0.005;
3592       fC[12][1] = 0.005;
3593       fC[12][2] = 1.0000;
3594       fC[12][3] = 0.3165;
3595       fC[12][4] = 0.3564;
3596
3597       fC[13][0] = 0.005;
3598       fC[13][1] = 0.005;
3599       fC[13][2] = 1.0000;
3600       fC[13][3] = 0.3526;
3601       fC[13][4] = 0.3011;
3602
3603       fC[14][0] = 0.005;
3604       fC[14][1] = 0.005;
3605       fC[14][2] = 1.0000;
3606       fC[14][3] = 0.3788;
3607       fC[14][4] = 0.3011;
3608
3609       fC[15][0] = 0.005;
3610       fC[15][1] = 0.005;
3611       fC[15][2] = 1.0000;
3612       fC[15][3] = 0.3788;
3613       fC[15][4] = 0.3011;
3614
3615       fC[16][0] = 0.005;
3616       fC[16][1] = 0.005;
3617       fC[16][2] = 1.0000;
3618       fC[16][3] = 0.3788;
3619       fC[16][4] = 0.3011;
3620
3621       fC[17][0] = 0.005;
3622       fC[17][1] = 0.005;
3623       fC[17][2] = 1.0000;
3624       fC[17][3] = 0.3788;
3625       fC[17][4] = 0.3011;
3626   }
3627
3628   for(Int_t i=18;i<fgkPIDptBin;i++){
3629     fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
3630     fC[i][0] = fC[17][0];
3631     fC[i][1] = fC[17][1];
3632     fC[i][2] = fC[17][2];
3633     fC[i][3] = fC[17][3];
3634     fC[i][4] = fC[17][4];
3635   }
3636 }
3637
3638 //---------------------------------------------------------------//
3639 Bool_t AliFlowTrackCuts::TPCTOFagree(const AliESDtrack *track)
3640 {
3641   //check pid agreement between TPC and TOF
3642   Bool_t status = kFALSE;
3643
3644   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
3645
3646
3647   Double_t exptimes[5];
3648   track->GetIntegratedTimes(exptimes);
3649
3650   Float_t dedx = track->GetTPCsignal();
3651
3652   Float_t p = track->P();
3653   Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
3654   Float_t tl = track->GetIntegratedLength();
3655
3656   Float_t betagammares =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
3657
3658   Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
3659
3660 //  printf("betagamma1 = %f\n",betagamma1);
3661
3662   if(betagamma1 < 0.1) betagamma1 = 0.1;
3663
3664   if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
3665   else betagamma1 = 100;
3666
3667   Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
3668 //  printf("betagamma2 = %f\n",betagamma2);
3669
3670   if(betagamma2 < 0.1) betagamma2 = 0.1;
3671
3672   if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
3673   else betagamma2 = 100;
3674
3675
3676   Double_t ptpc[3];
3677   track->GetInnerPxPyPz(ptpc);
3678   Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
3679
3680   for(Int_t i=0;i < 5;i++){
3681     Float_t resolutionTOF =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
3682     if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
3683       Float_t dedxExp = 0;
3684       if(i==0) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
3685       else if(i==1) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
3686       else if(i==2) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
3687       else if(i==3) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
3688       else if(i==4) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
3689
3690       Float_t resolutionTPC = 2;
3691       if(i==0) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
3692       else if(i==1) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
3693       else if(i==2) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
3694       else if(i==3) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
3695       else if(i==4) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
3696
3697       if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
3698         status = kTRUE;
3699       }
3700     }
3701   }
3702
3703   Float_t bb1 =  fESDpid.GetTPCResponse().Bethe(betagamma1);
3704   Float_t bb2 =  fESDpid.GetTPCResponse().Bethe(betagamma2);
3705   Float_t bbM =  fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
3706
3707
3708   //  status = kFALSE;
3709   // for nuclei
3710   Float_t resolutionTOFpr =   fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
3711   Float_t resolutionTPCpr =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
3712   if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3713      status = kTRUE;
3714   }
3715   else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3716      status = kTRUE;
3717   }
3718   else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
3719      status = kTRUE;
3720   }
3721
3722   return status;
3723 }
3724 // end part added by F. Noferini
3725 //-----------------------------------------------------------------------
3726
3727 //-----------------------------------------------------------------------
3728 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
3729 {
3730   //get the name of the particle id source
3731   switch (s)
3732   {
3733     case kTPCdedx:
3734       return "TPCdedx";
3735     case kTPCbayesian:
3736       return "TPCbayesian";
3737     case kTOFbeta:
3738       return "TOFbeta";
3739     case kTPCpid:
3740       return "TPCpid";
3741     case kTOFpid:
3742       return "TOFpid";
3743     case kTOFbayesian:
3744       return "TOFbayesianPID";
3745     case kTOFbetaSimple:
3746       return "TOFbetaSimple";
3747     case kTPCNuclei:
3748       return "TPCnuclei";
3749     default:
3750       return "NOPID";
3751   }
3752 }
3753
3754 //-----------------------------------------------------------------------
3755 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
3756 {
3757   //return the name of the selected parameter type
3758   switch (type)
3759   {
3760     case kMC:
3761       return "MC";
3762     case kGlobal:
3763       return "Global";
3764     case kTPCstandalone:
3765       return "TPCstandalone";
3766     case kSPDtracklet:
3767       return "SPDtracklets";
3768     case kPMD:
3769       return "PMD";
3770     case kV0:
3771       return "V0";
3772     case kMUON:       // XZhang 20120604
3773       return "MUON";  // XZhang 20120604
3774     default:
3775       return "unknown";
3776   }
3777 }
3778
3779 //-----------------------------------------------------------------------
3780 Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
3781 {
3782   //check PMD specific cuts
3783   //clean up from last iteration, and init label
3784   Int_t   det   = track->GetDetector();
3785   //Int_t   smn   = track->GetSmn();
3786   Float_t clsX  = track->GetClusterX();
3787   Float_t clsY  = track->GetClusterY();
3788   Float_t clsZ  = track->GetClusterZ();
3789   Float_t ncell = track->GetClusterCells();
3790   Float_t adc   = track->GetClusterADC();
3791
3792   fTrack = NULL;
3793   fMCparticle=NULL;
3794   fTrackLabel=-996;
3795
3796   fTrackEta = GetPmdEta(clsX,clsY,clsZ);
3797   fTrackPhi = GetPmdPhi(clsX,clsY);
3798   fTrackWeight = 1.0;
3799
3800   Bool_t pass=kTRUE;
3801   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
3802   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
3803   if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
3804   if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
3805   if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
3806
3807   return pass;
3808 }
3809
3810 //-----------------------------------------------------------------------
3811 Bool_t AliFlowTrackCuts::PassesV0cuts(Int_t id)
3812 {
3813   //check V0 cuts
3814   if (id<0) return kFALSE;
3815
3816   //clean up from last iter
3817   fTrack = NULL;
3818   fMCparticle=NULL;
3819   fTrackLabel=-995;
3820
3821   fTrackPhi = TMath::PiOver4()*(0.5+id%8);
3822
3823   if(id<32)
3824     fTrackEta = -3.45+0.5*(id/8);
3825   else
3826     fTrackEta = +4.8-0.6*((id/8)-4);
3827
3828   fTrackWeight = fEvent->GetVZEROEqMultiplicity(id);
3829
3830   if (fLinearizeVZEROresponse)
3831   {
3832     //this is only needed in pass1 of LHC10h
3833     Float_t multV0[fgkNumberOfV0tracks];
3834     Float_t dummy=0.0;
3835     AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multV0);
3836     fTrackWeight = multV0[id];
3837   }
3838
3839   Bool_t pass=kTRUE;
3840   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
3841   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
3842
3843   return pass;
3844 }
3845
3846 //-----------------------------------------------------------------------------
3847 Bool_t AliFlowTrackCuts::PassesMuonCuts(AliVParticle* vparticle)
3848 {
3849 // XZhang 20120604
3850   fTrack=NULL;
3851   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
3852   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
3853   else fMCparticle=NULL;
3854
3855   AliESDMuonTrack *esdTrack = dynamic_cast<AliESDMuonTrack*>(vparticle);
3856   AliAODTrack     *aodTrack = dynamic_cast<AliAODTrack*>(vparticle);
3857   if ((!esdTrack) && (!aodTrack)) return kFALSE;
3858   if (!fMuonTrackCuts->IsSelected(vparticle)) return kFALSE;
3859   HandleVParticle(vparticle);   if (!fTrack)  return kFALSE;
3860   return kTRUE;
3861 }
3862
3863 //----------------------------------------------------------------------------//
3864 Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
3865 {
3866   //get PMD "track" eta coordinate
3867   Float_t rpxpy, theta, eta;
3868   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
3869   theta  = TMath::ATan2(rpxpy,zPos);
3870   eta    = -TMath::Log(TMath::Tan(0.5*theta));
3871   return eta;
3872 }
3873
3874 //--------------------------------------------------------------------------//
3875 Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
3876 {
3877   //Get PMD "track" phi coordinate
3878   Float_t pybypx, phi = 0., phi1;
3879   if(xPos==0)
3880     {
3881       if(yPos>0) phi = 90.;
3882       if(yPos<0) phi = 270.;
3883     }
3884   if(xPos != 0)
3885     {
3886       pybypx = yPos/xPos;
3887       if(pybypx < 0) pybypx = - pybypx;
3888       phi1 = TMath::ATan(pybypx)*180./3.14159;
3889
3890       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
3891       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
3892       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
3893       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
3894
3895     }
3896   phi = phi*3.14159/180.;
3897   return   phi;
3898 }
3899
3900 //---------------------------------------------------------------//
3901 void AliFlowTrackCuts::Browse(TBrowser* b)
3902 {
3903   //some browsing capabilities
3904   if (fQA) b->Add(fQA);
3905 }
3906
3907 //---------------------------------------------------------------//
3908 Long64_t AliFlowTrackCuts::Merge(TCollection* list)
3909 {
3910   //merge
3911   if (!fQA || !list) return 0;
3912   if (list->IsEmpty()) return 0;
3913   AliFlowTrackCuts* obj=NULL;
3914   TList tmplist;
3915   TIter next(list);
3916   while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
3917   {
3918     if (obj==this) continue;
3919     tmplist.Add(obj->GetQA());
3920   }
3921   return fQA->Merge(&tmplist);
3922 }
3923
3924
3925