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