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 "TParticle.h"
  41 #include "TObjArray.h"
  42 #include "AliStack.h"
  43 #include "AliMCEvent.h"
  44 #include "AliESDEvent.h"
  45 #include "AliVParticle.h"
  46 #include "AliMCParticle.h"
  47 #include "AliESDtrack.h"
  48 #include "AliMultiplicity.h"
  49 #include "AliAODTrack.h"
  50 #include "AliFlowTrack.h"
  51 #include "AliFlowTrackCuts.h"
  52 #include "AliLog.h"
  53 #include "AliESDpid.h"
  54
  55 ClassImp(AliFlowTrackCuts)
  56
  57 //-----------------------------------------------------------------------
  58 AliFlowTrackCuts::AliFlowTrackCuts():
  59   AliFlowTrackSimpleCuts(),
  60   fAliESDtrackCuts(NULL),
  61   fQA(NULL),
  62   fCutMC(kFALSE),
  63   fCutMCprocessType(kFALSE),
  64   fMCprocessType(kPNoProcess),
  65   fCutMCPID(kFALSE),
  66   fMCPID(0),
  67   fIgnoreSignInPID(kFALSE),
  68   fCutMCisPrimary(kFALSE),
  69   fRequireTransportBitForPrimaries(kTRUE),
  70   fMCisPrimary(kFALSE),
  71   fRequireCharge(kFALSE),
  72   fFakesAreOK(kTRUE),
  73   fCutSPDtrackletDeltaPhi(kFALSE),
  74   fSPDtrackletDeltaPhiMax(FLT_MAX),
  75   fSPDtrackletDeltaPhiMin(-FLT_MAX),
  76   fIgnoreTPCzRange(kFALSE),
  77   fIgnoreTPCzRangeMax(FLT_MAX),
  78   fIgnoreTPCzRangeMin(-FLT_MAX),
  79   fCutChi2PerClusterTPC(kFALSE),
  80   fMaxChi2PerClusterTPC(FLT_MAX),
  81   fMinChi2PerClusterTPC(-FLT_MAX),
  82   fCutNClustersTPC(kFALSE),
  83   fNClustersTPCMax(INT_MAX),
  84   fNClustersTPCMin(INT_MIN),
  85   fCutNClustersITS(kFALSE),
  86   fNClustersITSMax(INT_MAX),
  87   fNClustersITSMin(INT_MIN),
  88   fUseAODFilterBit(kFALSE),
  89   fAODFilterBit(0),
  90   fCutDCAToVertexXY(kFALSE),
  91   fCutDCAToVertexZ(kFALSE),
  92   fCutMinimalTPCdedx(kFALSE),
  93   fMinimalTPCdedx(0.),
  94   fParamType(kGlobal),
  95   fParamMix(kPure),
  96   fTrack(NULL),
  97   fTrackPhi(0.),
  98   fTrackEta(0.),
  99   fTrackWeight(0.),
 100   fTrackLabel(INT_MIN),
 101   fMCevent(NULL),
 102   fMCparticle(NULL),
 103   fEvent(NULL),
 104   fTPCtrack(),
 105   fESDpid(),
 106   fPIDsource(kTOFpid),
 107   fTPCpidCuts(NULL),
 108   fTOFpidCuts(NULL),
 109   fParticleID(AliPID::kPion),
 110   fParticleProbability(.9)
 111 {
 112   //io constructor
 113   for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
 114   SetPriors(); //init arrays
 115 }
 116
 117 //-----------------------------------------------------------------------
 118 AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
 119   AliFlowTrackSimpleCuts(),
 120   fAliESDtrackCuts(new AliESDtrackCuts()),
 121   fQA(NULL),
 122   fCutMC(kFALSE),
 123   fCutMCprocessType(kFALSE),
 124   fMCprocessType(kPNoProcess),
 125   fCutMCPID(kFALSE),
 126   fMCPID(0),
 127   fIgnoreSignInPID(kFALSE),
 128   fCutMCisPrimary(kFALSE),
 129   fRequireTransportBitForPrimaries(kTRUE),
 130   fMCisPrimary(kFALSE),
 131   fRequireCharge(kFALSE),
 132   fFakesAreOK(kTRUE),
 133   fCutSPDtrackletDeltaPhi(kFALSE),
 134   fSPDtrackletDeltaPhiMax(FLT_MAX),
 135   fSPDtrackletDeltaPhiMin(-FLT_MAX),
 136   fIgnoreTPCzRange(kFALSE),
 137   fIgnoreTPCzRangeMax(FLT_MAX),
 138   fIgnoreTPCzRangeMin(-FLT_MAX),
 139   fCutChi2PerClusterTPC(kFALSE),
 140   fMaxChi2PerClusterTPC(FLT_MAX),
 141   fMinChi2PerClusterTPC(-FLT_MAX),
 142   fCutNClustersTPC(kFALSE),
 143   fNClustersTPCMax(INT_MAX),
 144   fNClustersTPCMin(INT_MIN),
 145   fCutNClustersITS(kFALSE),
 146   fNClustersITSMax(INT_MAX),
 147   fNClustersITSMin(INT_MIN),
 148   fUseAODFilterBit(kFALSE),
 149   fAODFilterBit(0),
 150   fCutDCAToVertexXY(kFALSE),
 151   fCutDCAToVertexZ(kFALSE),
 152   fCutMinimalTPCdedx(kFALSE),
 153   fMinimalTPCdedx(0.),
 154   fParamType(kGlobal),
 155   fParamMix(kPure),
 156   fTrack(NULL),
 157   fTrackPhi(0.),
 158   fTrackEta(0.),
 159   fTrackWeight(0.),
 160   fTrackLabel(INT_MIN),
 161   fMCevent(NULL),
 162   fMCparticle(NULL),
 163   fEvent(NULL),
 164   fTPCtrack(),
 165   fESDpid(),
 166   fPIDsource(kTOFpid),
 167   fTPCpidCuts(NULL),
 168   fTOFpidCuts(NULL),
 169   fParticleID(AliPID::kPion),
 170   fParticleProbability(.9)
 171 {
 172   //constructor
 173   SetName(name);
 174   SetTitle("AliFlowTrackCuts");
 175   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
 176                                                     2.63394e+01,
 177                                                     5.04114e-11,
 178                                                     2.12543e+00,
 179                                                     4.88663e+00 );
 180   for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
 181   SetPriors(); //init arrays
 182 }
 183
 184 //-----------------------------------------------------------------------
 185 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
 186   AliFlowTrackSimpleCuts(that),
 187   fAliESDtrackCuts(NULL),
 188   fQA(NULL),
 189   fCutMC(that.fCutMC),
 190   fCutMCprocessType(that.fCutMCprocessType),
 191   fMCprocessType(that.fMCprocessType),
 192   fCutMCPID(that.fCutMCPID),
 193   fMCPID(that.fMCPID),
 194   fIgnoreSignInPID(that.fIgnoreSignInPID),
 195   fCutMCisPrimary(that.fCutMCisPrimary),
 196   fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
 197   fMCisPrimary(that.fMCisPrimary),
 198   fRequireCharge(that.fRequireCharge),
 199   fFakesAreOK(that.fFakesAreOK),
 200   fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
 201   fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
 202   fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
 203   fIgnoreTPCzRange(that.fIgnoreTPCzRange),
 204   fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
 205   fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
 206   fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
 207   fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
 208   fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
 209   fCutNClustersTPC(that.fCutNClustersTPC),
 210   fNClustersTPCMax(that.fNClustersTPCMax),
 211   fNClustersTPCMin(that.fNClustersTPCMin),
 212   fCutNClustersITS(that.fCutNClustersITS),
 213   fNClustersITSMax(that.fNClustersITSMax),
 214   fNClustersITSMin(that.fNClustersITSMin),
 215   fUseAODFilterBit(that.fUseAODFilterBit),
 216   fAODFilterBit(that.fAODFilterBit),
 217   fCutDCAToVertexXY(that.fCutDCAToVertexXY),
 218   fCutDCAToVertexZ(that.fCutDCAToVertexZ),
 219   fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
 220   fMinimalTPCdedx(that.fMinimalTPCdedx),
 221   fParamType(that.fParamType),
 222   fParamMix(that.fParamMix),
 223   fTrack(NULL),
 224   fTrackPhi(0.),
 225   fTrackEta(0.),
 226   fTrackWeight(0.),
 227   fTrackLabel(INT_MIN),
 228   fMCevent(NULL),
 229   fMCparticle(NULL),
 230   fEvent(NULL),
 231   fTPCtrack(),
 232   fESDpid(that.fESDpid),
 233   fPIDsource(that.fPIDsource),
 234   fTPCpidCuts(NULL),
 235   fTOFpidCuts(NULL),
 236   fParticleID(that.fParticleID),
 237   fParticleProbability(that.fParticleProbability)
 238 {
 239   //copy constructor
 240   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 241   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 242   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
 243   memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
 244   SetPriors(); //init arrays
 245 }
 246
 247 //-----------------------------------------------------------------------
 248 AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
 249 {
 250   //assignment
 251   if (this==&that) return *this;
 252
 253   AliFlowTrackSimpleCuts::operator=(that);
 254   if (that.fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
 255   fQA=NULL;
 256   fCutMC=that.fCutMC;
 257   fCutMCprocessType=that.fCutMCprocessType;
 258   fMCprocessType=that.fMCprocessType;
 259   fCutMCPID=that.fCutMCPID;
 260   fMCPID=that.fMCPID;
 261   fIgnoreSignInPID=that.fIgnoreSignInPID,
 262   fCutMCisPrimary=that.fCutMCisPrimary;
 263   fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
 264   fMCisPrimary=that.fMCisPrimary;
 265   fRequireCharge=that.fRequireCharge;
 266   fFakesAreOK=that.fFakesAreOK;
 267   fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
 268   fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
 269   fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
 270   fIgnoreTPCzRange=that.fIgnoreTPCzRange;
 271   fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
 272   fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
 273   fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
 274   fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
 275   fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
 276   fCutNClustersTPC=that.fCutNClustersTPC;
 277   fNClustersTPCMax=that.fNClustersTPCMax;
 278   fNClustersTPCMin=that.fNClustersTPCMin;
 279   fCutNClustersITS=that.fCutNClustersITS;
 280   fNClustersITSMax=that.fNClustersITSMax;
 281   fNClustersITSMin=that.fNClustersITSMin;
 282   fUseAODFilterBit=that.fUseAODFilterBit;
 283   fAODFilterBit=that.fAODFilterBit;
 284   fCutDCAToVertexXY=that.fCutDCAToVertexXY;
 285   fCutDCAToVertexZ=that.fCutDCAToVertexZ;
 286   fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
 287   fMinimalTPCdedx=that.fMinimalTPCdedx;
 288   fParamType=that.fParamType;
 289   fParamMix=that.fParamMix;
 290
 291   fTrack=NULL;
 292   fTrackPhi=0.;
 293   fTrackPhi=0.;
 294   fTrackWeight=0.;
 295   fTrackLabel=INT_MIN;
 296   fMCevent=NULL;
 297   fMCparticle=NULL;
 298   fEvent=NULL;
 299
 300   fESDpid = that.fESDpid;
 301   fPIDsource = that.fPIDsource;
 302
 303   delete fTPCpidCuts;
 304   delete fTOFpidCuts;
 305   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 306   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 307
 308   fParticleID=that.fParticleID;
 309   fParticleProbability=that.fParticleProbability;
 310   memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
 311
 312   return *this;
 313 }
 314
 315 //-----------------------------------------------------------------------
 316 AliFlowTrackCuts::~AliFlowTrackCuts()
 317 {
 318   //dtor
 319   delete fAliESDtrackCuts;
 320   delete fTPCpidCuts;
 321   delete fTOFpidCuts;
 322 }
 323
 324 //-----------------------------------------------------------------------
 325 void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
 326 {
 327   //set the event
 328   Clear();
 329   fEvent=event;
 330   fMCevent=mcEvent;
 331
 332   //do the magic for ESD
 333   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
 334   if (fCutPID && myESD)
 335   {
 336     //TODO: maybe call it only for the TOF options?
 337     // Added by F. Noferini for TOF PID
 338     fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
 339     fESDpid.MakePID(myESD,kFALSE);
 340     // End F. Noferini added part
 341   }
 342
 343   //TODO: AOD
 344 }
 345
 346 //-----------------------------------------------------------------------
 347 void AliFlowTrackCuts::SetCutMC( Bool_t b )
 348 {
 349   //will we be cutting on MC information?
 350   fCutMC=b;
 351
 352   //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
 353   if (fCutMC)
 354   {
 355     fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
 356                                                        1.75295e+01,
 357                                                        3.40030e-09,
 358                                                        1.96178e+00,
 359                                                        3.91720e+00);
 360   }
 361 }
 362
 363 //-----------------------------------------------------------------------
 364 Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
 365 {
 366   //check cuts
 367   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 368   if (vparticle) return PassesCuts(vparticle);
 369   AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
 370   if (flowtrack) return PassesCuts(flowtrack);
 371   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 372   if (tracklets) return PassesCuts(tracklets,id);
 373   return kFALSE;  //default when passed wrong type of object
 374 }
 375
 376 //-----------------------------------------------------------------------
 377 Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
 378 {
 379   //check cuts
 380   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 381   if (vparticle)
 382   {
 383     return PassesMCcuts(fMCevent,vparticle->GetLabel());
 384   }
 385   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 386   if (tracklets)
 387   {
 388     Int_t label0 = tracklets->GetLabel(id,0);
 389     Int_t label1 = tracklets->GetLabel(id,1);
 390     Int_t label = (label0==label1)?tracklets->GetLabel(id,1):-666;
 391     return PassesMCcuts(fMCevent,label);
 392   }
 393   return kFALSE;  //default when passed wrong type of object
 394 }
 395
 396 //-----------------------------------------------------------------------
 397 Bool_t AliFlowTrackCuts::PassesCuts(AliFlowTrackSimple* track)
 398 {
 399   //check cuts on a flowtracksimple
 400
 401   //clean up from last iteration
 402   fTrack = NULL;
 403   return AliFlowTrackSimpleCuts::PassesCuts(track);
 404 }
 405
 406 //-----------------------------------------------------------------------
 407 Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
 408 {
 409   //check cuts on a tracklets
 410
 411   //clean up from last iteration, and init label
 412   fTrack = NULL;
 413   fMCparticle=NULL;
 414   fTrackLabel=-1;
 415
 416   fTrackPhi = tracklet->GetPhi(id);
 417   fTrackEta = tracklet->GetEta(id);
 418   fTrackWeight = 1.0;
 419   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
 420   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
 421
 422   //check MC info if available
 423   //if the 2 clusters have different label track cannot be good
 424   //and should therefore not pass the mc cuts
 425   Int_t label0 = tracklet->GetLabel(id,0);
 426   Int_t label1 = tracklet->GetLabel(id,1);
 427   //if possible get label and mcparticle
 428   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
 429   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
 430   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 431   //check MC cuts
 432   if (fCutMC && !PassesMCcuts()) return kFALSE;
 433   return kTRUE;
 434 }
 435
 436 //-----------------------------------------------------------------------
 437 Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
 438 {
 439   //check the MC info
 440   if (!mcEvent) return kFALSE;
 441   if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 442   AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
 443   if (!mcparticle) {AliError("no MC track"); return kFALSE;}
 444
 445   if (fCutMCisPrimary)
 446   {
 447     if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
 448   }
 449   if (fCutMCPID)
 450   {
 451     Int_t pdgCode = mcparticle->PdgCode();
 452     if (fIgnoreSignInPID)
 453     {
 454       if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
 455     }
 456     else
 457     {
 458       if (fMCPID != pdgCode) return kFALSE;
 459     }
 460   }
 461   if ( fCutMCprocessType )
 462   {
 463     TParticle* particle = mcparticle->Particle();
 464     Int_t processID = particle->GetUniqueID();
 465     if (processID != fMCprocessType ) return kFALSE;
 466   }
 467   return kTRUE;
 468 }
 469 //-----------------------------------------------------------------------
 470 Bool_t AliFlowTrackCuts::PassesMCcuts()
 471 {
 472   if (!fMCevent) return kFALSE;
 473   if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 474   fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 475   return PassesMCcuts(fMCevent,fTrackLabel);
 476 }
 477
 478 //-----------------------------------------------------------------------
 479 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
 480 {
 481   //check cuts for an ESD vparticle
 482
 483   ////////////////////////////////////////////////////////////////
 484   //  start by preparing the track parameters to cut on //////////
 485   ////////////////////////////////////////////////////////////////
 486   //clean up from last iteration
 487   fTrack=NULL;
 488
 489   //get the label and the mc particle
 490   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
 491   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 492   else fMCparticle=NULL;
 493
 494   Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
 495   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
 496   AliAODTrack* aodTrack = NULL;
 497   if (esdTrack)
 498     //for an ESD track we do some magic sometimes like constructing TPC only parameters
 499     //or doing some hybrid, handle that here
 500     HandleESDtrack(esdTrack);
 501   else
 502   {
 503     HandleVParticle(vparticle);
 504     //now check if produced particle is MC
 505     isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
 506     aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
 507   }
 508   ////////////////////////////////////////////////////////////////
 509   ////////////////////////////////////////////////////////////////
 510
 511   if (!fTrack) return kFALSE;
 512   //because it may be different from global, not needed for aodTrack because we dont do anything funky there
 513   if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
 514
 515   Bool_t pass=kTRUE;
 516   //check the common cuts for the current particle fTrack (MC,AOD,ESD)
 517   Double_t pt = fTrack->Pt();
 518   if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
 519   if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
 520   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
 521   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
 522   if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
 523   if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
 524   if (fCutCharge && isMCparticle)
 525   {
 526     //in case of an MC particle the charge is stored in units of 1/3|e|
 527     Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
 528     if (charge!=fCharge) pass=kFALSE;
 529   }
 530   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
 531
 532   //when additionally MC info is required
 533   if (fCutMC && !PassesMCcuts()) pass=kFALSE;
 534
 535   //the case of ESD or AOD
 536   if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
 537   if (aodTrack) { if (!PassesAODcuts(aodTrack)) { pass=kFALSE; } }
 538
 539   //true by default, if we didn't set any cuts
 540   return pass;
 541 }
 542
 543 //_______________________________________________________________________
 544 Bool_t AliFlowTrackCuts::PassesAODcuts(AliAODTrack* track)
 545 {
 546   Bool_t pass = kTRUE;
 547
 548   if (fCutNClustersTPC)
 549   {
 550     Int_t ntpccls = track->GetTPCNcls();
 551     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
 552   }
 553
 554   if (fCutNClustersITS)
 555   {
 556     Int_t nitscls = track->GetITSNcls();
 557     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
 558   }
 559
 560    if (fCutChi2PerClusterTPC)
 561   {
 562     Double_t chi2tpc = track->Chi2perNDF();
 563     if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
 564   }
 565
 566   if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
 567   if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
 568
 569   if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
 570
 571   if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
 572
 573
 574   return pass;
 575 }
 576
 577 //_______________________________________________________________________
 578 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
 579 {
 580   Bool_t pass=kTRUE;
 581   if (fIgnoreTPCzRange)
 582   {
 583     const AliExternalTrackParam* pin = track->GetOuterParam();
 584     const AliExternalTrackParam* pout = track->GetInnerParam();
 585     if (pin&&pout)
 586     {
 587       Double_t zin = pin->GetZ();
 588       Double_t zout = pout->GetZ();
 589       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
 590       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
 591       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
 592     }
 593   }
 594
 595   Int_t ntpccls = ( fParamType==kESD_TPConly )?
 596                     track->GetTPCNclsIter1():track->GetTPCNcls();
 597   if (fCutChi2PerClusterTPC)
 598   {
 599     Float_t tpcchi2 = (fParamType==kESD_TPConly)?
 600                        track->GetTPCchi2Iter1():track->GetTPCchi2();
 601     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
 602     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
 603       pass=kFALSE;
 604   }
 605
 606   if (fCutMinimalTPCdedx)
 607   {
 608     if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
 609   }
 610
 611   if (fCutNClustersTPC)
 612   {
 613     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
 614   }
 615
 616   Int_t nitscls = track->GetNcls(0);
 617   if (fCutNClustersITS)
 618   {
 619     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
 620   }
 621
 622   if (fCutPID)
 623   {
 624     switch (fPIDsource)
 625     {
 626       case kTPCpid:
 627         if (!PassesTPCpidCut(track)) pass=kFALSE;
 628         break;
 629       case kTPCdedx:
 630         if (!PassesTPCdedxCut(track)) pass=kFALSE;
 631         break;
 632       case kTOFpid:
 633         if (!PassesTOFpidCut(track)) pass=kFALSE;
 634         break;
 635       case kTOFbeta:
 636         if (!PassesTOFbetaCut(track)) pass=kFALSE;
 637         break;
 638             // part added by F. Noferini
 639       case kTOFbayesian:
 640               if (!PassesTOFbayesianCut(track)) pass=kFALSE;
 641               break;
 642             // end part added by F. Noferini
 643       default:
 644         printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
 645         pass=kFALSE;
 646         break;
 647     }
 648   }
 649
 650   //some stuff is still handled by AliESDtrackCuts class - delegate
 651   if (fAliESDtrackCuts)
 652   {
 653     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
 654   }
 655
 656   return pass;
 657 }
 658
 659 //-----------------------------------------------------------------------
 660 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
 661 {
 662   //handle the general case
 663   switch (fParamType)
 664   {
 665     default:
 666       fTrack = track;
 667       break;
 668   }
 669 }
 670
 671 //-----------------------------------------------------------------------
 672 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
 673 {
 674   //handle esd track
 675   switch (fParamType)
 676   {
 677     case kGlobal:
 678       fTrack = track;
 679       break;
 680     case kESD_TPConly:
 681       if (!track->FillTPCOnlyTrack(fTPCtrack))
 682       {
 683         fTrack=NULL;
 684         fMCparticle=NULL;
 685         fTrackLabel=-1;
 686         return;
 687       }
 688       fTrack = &fTPCtrack;
 689       //recalculate the label and mc particle, they may differ as TPClabel != global label
 690       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
 691       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 692       else fMCparticle=NULL;
 693       break;
 694     default:
 695       fTrack = track;
 696       break;
 697   }
 698 }
 699
 700 //-----------------------------------------------------------------------
 701 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
 702 {
 703   //calculate the number of track in given event.
 704   //if argument is NULL(default) take the event attached
 705   //by SetEvent()
 706   Int_t multiplicity = 0;
 707   if (!event)
 708   {
 709     for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
 710     {
 711       if (IsSelected(GetInputObject(i))) multiplicity++;
 712     }
 713   }
 714   else
 715   {
 716     for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
 717     {
 718       if (IsSelected(event->GetTrack(i))) multiplicity++;
 719     }
 720   }
 721   return multiplicity;
 722 }
 723
 724 //-----------------------------------------------------------------------
 725 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
 726 {
 727   //get standard cuts
 728   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global flow cuts");
 729   cuts->SetParamType(kGlobal);
 730   cuts->SetPtRange(0.2,5.);
 731   cuts->SetEtaRange(-0.8,0.8);
 732   cuts->SetMinNClustersTPC(70);
 733   cuts->SetMinChi2PerClusterTPC(0.1);
 734   cuts->SetMaxChi2PerClusterTPC(4.0);
 735   cuts->SetMinNClustersITS(2);
 736   cuts->SetRequireITSRefit(kTRUE);
 737   cuts->SetRequireTPCRefit(kTRUE);
 738   cuts->SetMaxDCAToVertexXY(0.3);
 739   cuts->SetMaxDCAToVertexZ(0.3);
 740   cuts->SetAcceptKinkDaughters(kFALSE);
 741   cuts->SetMinimalTPCdedx(10.);
 742
 743   return cuts;
 744 }
 745
 746 //-----------------------------------------------------------------------
 747 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts2010()
 748 {
 749   //get standard cuts
 750   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly flow cuts");
 751   cuts->SetParamType(kESD_TPConly);
 752   cuts->SetPtRange(0.2,5.);
 753   cuts->SetEtaRange(-0.8,0.8);
 754   cuts->SetMinNClustersTPC(70);
 755   cuts->SetMinChi2PerClusterTPC(0.2);
 756   cuts->SetMaxChi2PerClusterTPC(4.0);
 757   cuts->SetMaxDCAToVertexXY(3.0);
 758   cuts->SetMaxDCAToVertexZ(3.0);
 759   cuts->SetDCAToVertex2D(kTRUE);
 760   cuts->SetAcceptKinkDaughters(kFALSE);
 761   cuts->SetMinimalTPCdedx(10.);
 762
 763   return cuts;
 764 }
 765
 766 //-----------------------------------------------------------------------
 767 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
 768 {
 769   //get standard cuts
 770   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly flow cuts");
 771   cuts->SetParamType(kESD_TPConly);
 772   cuts->SetPtRange(0.2,5.);
 773   cuts->SetEtaRange(-0.8,0.8);
 774   cuts->SetMinNClustersTPC(70);
 775   cuts->SetMinChi2PerClusterTPC(0.2);
 776   cuts->SetMaxChi2PerClusterTPC(4.0);
 777   cuts->SetMaxDCAToVertexXY(3.0);
 778   cuts->SetMaxDCAToVertexZ(3.0);
 779   cuts->SetDCAToVertex2D(kTRUE);
 780   cuts->SetAcceptKinkDaughters(kFALSE);
 781   cuts->SetMinimalTPCdedx(10.);
 782
 783   return cuts;
 784 }
 785
 786 //-----------------------------------------------------------------------
 787 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
 788 {
 789   //get standard cuts
 790   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
 791   delete cuts->fAliESDtrackCuts;
 792   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
 793   cuts->SetParamType(kGlobal);
 794   return cuts;
 795 }
 796
 797 //-----------------------------------------------------------------------
 798 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
 799 {
 800   //get a flow track constructed from whatever we applied cuts on
 801   //caller is resposible for deletion
 802   //if construction fails return NULL
 803   AliFlowTrack* flowtrack=NULL;
 804   TParticle *tmpTParticle=NULL;
 805   AliMCParticle* tmpAliMCParticle=NULL;
 806   if (fParamType==kESD_SPDtracklet)
 807   {
 808     switch (fParamMix)
 809     {
 810       case kPure:
 811         flowtrack = new AliFlowTrack();
 812         flowtrack->SetPhi(fTrackPhi);
 813         flowtrack->SetEta(fTrackEta);
 814         break;
 815       case kTrackWithMCkine:
 816         if (!fMCparticle) return NULL;
 817         flowtrack = new AliFlowTrack();
 818         flowtrack->SetPhi( fMCparticle->Phi() );
 819         flowtrack->SetEta( fMCparticle->Eta() );
 820         flowtrack->SetPt( fMCparticle->Pt() );
 821         break;
 822       case kTrackWithMCpt:
 823         if (!fMCparticle) return NULL;
 824         flowtrack = new AliFlowTrack();
 825         flowtrack->SetPhi(fTrackPhi);
 826         flowtrack->SetEta(fTrackEta);
 827         flowtrack->SetPt(fMCparticle->Pt());
 828         break;
 829       case kTrackWithPtFromFirstMother:
 830         if (!fMCparticle) return NULL;
 831         flowtrack = new AliFlowTrack();
 832         flowtrack->SetPhi(fTrackPhi);
 833         flowtrack->SetEta(fTrackEta);
 834         tmpTParticle = fMCparticle->Particle();
 835         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
 836         flowtrack->SetPt(tmpAliMCParticle->Pt());
 837         break;
 838       default:
 839         flowtrack = new AliFlowTrack();
 840         flowtrack->SetPhi(fTrackPhi);
 841         flowtrack->SetEta(fTrackEta);
 842         break;
 843     }
 844     flowtrack->SetSource(AliFlowTrack::kFromTracklet);
 845   }
 846   else
 847   {
 848     if (!fTrack) return NULL;
 849     switch(fParamMix)
 850     {
 851       case kPure:
 852         flowtrack = new AliFlowTrack(fTrack);
 853         break;
 854       case kTrackWithMCkine:
 855         flowtrack = new AliFlowTrack(fMCparticle);
 856         break;
 857       case kTrackWithMCPID:
 858         flowtrack = new AliFlowTrack(fTrack);
 859         //flowtrack->setPID(...) from mc, when implemented
 860         break;
 861       case kTrackWithMCpt:
 862         if (!fMCparticle) return NULL;
 863         flowtrack = new AliFlowTrack(fTrack);
 864         flowtrack->SetPt(fMCparticle->Pt());
 865         break;
 866       case kTrackWithPtFromFirstMother:
 867         if (!fMCparticle) return NULL;
 868         flowtrack = new AliFlowTrack(fTrack);
 869         tmpTParticle = fMCparticle->Particle();
 870         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
 871         flowtrack->SetPt(tmpAliMCParticle->Pt());
 872         break;
 873       default:
 874         flowtrack = new AliFlowTrack(fTrack);
 875         break;
 876     }
 877     if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
 878     else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
 879     else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
 880     else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
 881   }
 882
 883   return flowtrack;
 884 }
 885
 886 //-----------------------------------------------------------------------
 887 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
 888 {
 889   //check if current particle is a physical primary
 890   if (!fMCevent) return kFALSE;
 891   if (fTrackLabel<0) return kFALSE;
 892   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
 893 }
 894
 895 //-----------------------------------------------------------------------
 896 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
 897 {
 898   //check if current particle is a physical primary
 899   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
 900   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
 901   if (!track) return kFALSE;
 902   TParticle* particle = track->Particle();
 903   Bool_t transported = particle->TestBit(kTransportBit);
 904   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
 905         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
 906   return (physprim && (transported || !requiretransported));
 907 }
 908
 909 //-----------------------------------------------------------------------
 910 void AliFlowTrackCuts::DefineHistograms()
 911 {
 912   //define qa histograms
 913 }
 914
 915 //-----------------------------------------------------------------------
 916 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
 917 {
 918   //get the number of tracks in the input event according source
 919   //selection (ESD tracks, tracklets, MC particles etc.)
 920   AliESDEvent* esd=NULL;
 921   switch (fParamType)
 922   {
 923     case kESD_SPDtracklet:
 924       esd = dynamic_cast<AliESDEvent*>(fEvent);
 925       if (!esd) return 0;
 926       return esd->GetMultiplicity()->GetNumberOfTracklets();
 927     case kMC:
 928       if (!fMCevent) return 0;
 929       return fMCevent->GetNumberOfTracks();
 930     default:
 931       if (!fEvent) return 0;
 932       return fEvent->GetNumberOfTracks();
 933   }
 934   return 0;
 935 }
 936
 937 //-----------------------------------------------------------------------
 938 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
 939 {
 940   //get the input object according the data source selection:
 941   //(esd tracks, traclets, mc particles,etc...)
 942   AliESDEvent* esd=NULL;
 943   switch (fParamType)
 944   {
 945     case kESD_SPDtracklet:
 946       esd = dynamic_cast<AliESDEvent*>(fEvent);
 947       if (!esd) return NULL;
 948       return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
 949     case kMC:
 950       if (!fMCevent) return NULL;
 951       return fMCevent->GetTrack(i);
 952     default:
 953       if (!fEvent) return NULL;
 954       return fEvent->GetTrack(i);
 955   }
 956 }
 957
 958 //-----------------------------------------------------------------------
 959 void AliFlowTrackCuts::Clear(Option_t*)
 960 {
 961   //clean up
 962   fTrack=NULL;
 963   fMCevent=NULL;
 964   fMCparticle=NULL;
 965   fTrackLabel=0;
 966   fTrackWeight=0.0;
 967   fTrackEta=0.0;
 968   fTrackPhi=0.0;
 969 }
 970
 971 //-----------------------------------------------------------------------
 972 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(AliESDtrack* track )
 973 {
 974   //check if passes PID cut using timing in TOF
 975   Bool_t goodtrack = (track) &&
 976                      (track->GetStatus() & AliESDtrack::kTOFpid) &&
 977                      (track->GetTOFsignal() > 12000) &&
 978                      (track->GetTOFsignal() < 100000) &&
 979                      (track->GetIntegratedLength() > 365) &&
 980                     !(track->GetStatus() & AliESDtrack::kTOFmismatch);
 981
 982   if (!goodtrack) return kFALSE;
 983
 984   const Float_t c = 2.99792457999999984e-02;
 985   Float_t p = track->GetP();
 986   Float_t L = track->GetIntegratedLength();
 987   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
 988   Float_t timeTOF = track->GetTOFsignal()- trackT0;
 989   Float_t beta = L/timeTOF/c;
 990   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
 991   track->GetIntegratedTimes(integratedTimes);
 992
 993   //construct the pid index because it's not AliPID::EParticleType
 994   Int_t pid = 0;
 995   switch (fParticleID)
 996   {
 997     case AliPID::kPion:
 998       pid=2;
 999       break;
1000     case AliPID::kKaon:
1001       pid=3;
1002       break;
1003     case AliPID::kProton:
1004       pid=4;
1005       break;
1006     default:
1007       return kFALSE;
1008   }
1009
1010   //signal to cut on
1011   Float_t s = beta-L/integratedTimes[pid]/c;
1012
1013   Float_t* arr = fTOFpidCuts->GetMatrixArray();
1014   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
1015   if (col<0) return kFALSE;
1016   Float_t min = (*fTOFpidCuts)(1,col);
1017   Float_t max = (*fTOFpidCuts)(2,col);
1018
1019   //printf("--------------TOF beta cut %s\n",(s>min && s<max)?"PASS":"FAIL");
1020   return (s>min && s<max);
1021 }
1022
1023 //-----------------------------------------------------------------------
1024 Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* track)
1025 {
1026   //check if passes PID cut using default TOF pid
1027   Double_t pidTOF[AliPID::kSPECIES];
1028   track->GetTOFpid(pidTOF);
1029   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
1030   return kFALSE;
1031 }
1032
1033 //-----------------------------------------------------------------------
1034 Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
1035 {
1036   //check if passes PID cut using default TPC pid
1037   Double_t pidTPC[AliPID::kSPECIES];
1038   track->GetTPCpid(pidTPC);
1039   Double_t probablity = 0.;
1040   switch (fParticleID)
1041   {
1042     case AliPID::kPion:
1043       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
1044       break;
1045     default:
1046       probablity = pidTPC[fParticleID];
1047   }
1048   if (probablity >= fParticleProbability) return kTRUE;
1049   return kFALSE;
1050 }
1051
1052 //-----------------------------------------------------------------------
1053 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(AliESDtrack* track)
1054 {
1055   //check if passes PID cut using dedx signal in the TPC
1056   if (!fTPCpidCuts)
1057   {
1058     printf("no TPCpidCuts\n");
1059     return kFALSE;
1060   }
1061
1062   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
1063   if (!tpcparam) return kFALSE;
1064   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fParticleID);
1065   Float_t sigTPC = track->GetTPCsignal();
1066   Float_t s = (sigTPC-sigExp)/sigExp;
1067   Double_t pt = track->Pt();
1068
1069   Float_t* arr = fTPCpidCuts->GetMatrixArray();
1070   Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
1071   if (col<0) return kFALSE;
1072   Float_t min = (*fTPCpidCuts)(1,col);
1073   Float_t max = (*fTPCpidCuts)(2,col);
1074
1075   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
1076   return (s>min && s<max);
1077 }
1078
1079 //-----------------------------------------------------------------------
1080 void AliFlowTrackCuts::InitPIDcuts()
1081 {
1082   //init matrices with PID cuts
1083   TMatrixF* t = NULL;
1084   if (!fTPCpidCuts)
1085   {
1086     if (fParticleID==AliPID::kPion)
1087     {
1088       t = new TMatrixF(3,10);
1089       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.2;
1090       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.2;
1091       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.25;
1092       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.25;
1093       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
1094       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
1095       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
1096       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
1097       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
1098       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  =     0;
1099     }
1100     else
1101     if (fParticleID==AliPID::kKaon)
1102     {
1103       t = new TMatrixF(3,7);
1104       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.4;
1105       (*t)(0,1)  = 0.25;  (*t)(1,1)  =-0.15;  (*t)(2,1)  = 0.4;
1106       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.1;  (*t)(2,2)  = 0.4;
1107       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.1;  (*t)(2,3)  = 0.4;
1108       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.6;
1109       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.6;
1110       (*t)(0,6)  = 0.50;  (*t)(1,6)  =    0;  (*t)(2,6)  =   0;
1111     }
1112     else
1113     if (fParticleID==AliPID::kProton)
1114     {
1115       t = new TMatrixF(3,16);
1116       (*t)(0,0)  = 0.20;  (*t)(1,0)  =     0;  (*t)(2,0)  =    0;
1117       (*t)(0,1)  = 0.25;  (*t)(1,1)  =  -0.2;  (*t)(2,1)  =  0.3;
1118       (*t)(0,2)  = 0.30;  (*t)(1,2)  =  -0.2;  (*t)(2,2)  =  0.6;
1119       (*t)(0,3)  = 0.35;  (*t)(1,3)  =  -0.2;  (*t)(2,3)  =  0.6;
1120       (*t)(0,4)  = 0.40;  (*t)(1,4)  =  -0.2;  (*t)(2,4)  =  0.6;
1121       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.15;  (*t)(2,5)  =  0.6;
1122       (*t)(0,6)  = 0.50;  (*t)(1,6)  =  -0.1;  (*t)(2,6)  =  0.6;
1123       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =  0.6;
1124       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.05;  (*t)(2,8)  = 0.45;
1125       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.05;  (*t)(2,9)  = 0.45;
1126       (*t)(0,10) = 0.70;  (*t)(1,10) = -0.05;  (*t)(2,10) = 0.45;
1127       (*t)(0,11) = 0.75;  (*t)(1,11) = -0.05;  (*t)(2,11) = 0.45;
1128       (*t)(0,12) = 0.80;  (*t)(1,12) =     0;  (*t)(2,12) = 0.45;
1129       (*t)(0,13) = 0.85;  (*t)(1,13) =     0;  (*t)(2,13) = 0.45;
1130       (*t)(0,14) = 0.90;  (*t)(1,14) =     0;  (*t)(2,14) = 0.45;
1131       (*t)(0,15) = 0.95;  (*t)(1,15) =     0;  (*t)(2,15) =    0;
1132     }
1133     fTPCpidCuts=t;
1134   }
1135   t = NULL;
1136   if (!fTOFpidCuts)
1137   {
1138     if (fParticleID==AliPID::kPion)
1139     {
1140       //TOF pions, 0.9 purity
1141       t = new TMatrixF(3,61);
1142       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1143       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1144       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1145       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1146       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1147       (*t)(0,5)  = 0.250;  (*t)(2,5)  = -0.046;  (*t)(2,5)  =   0.046;
1148       (*t)(0,6)  = 0.300;  (*t)(2,6)  = -0.038;  (*t)(2,6)  =   0.038;
1149       (*t)(0,7)  = 0.350;  (*t)(2,7)  = -0.034;  (*t)(2,7)  =   0.034;
1150       (*t)(0,8)  = 0.400;  (*t)(2,8)  = -0.032;  (*t)(2,8)  =   0.032;
1151       (*t)(0,9)  = 0.450;  (*t)(2,9)  = -0.030;  (*t)(2,9)  =   0.030;
1152       (*t)(0,10)  = 0.500;  (*t)(2,10)  = -0.030;  (*t)(2,10)  =   0.030;
1153       (*t)(0,11)  = 0.550;  (*t)(2,11)  = -0.030;  (*t)(2,11)  =   0.030;
1154       (*t)(0,12)  = 0.600;  (*t)(2,12)  = -0.030;  (*t)(2,12)  =   0.030;
1155       (*t)(0,13)  = 0.650;  (*t)(2,13)  = -0.030;  (*t)(2,13)  =   0.030;
1156       (*t)(0,14)  = 0.700;  (*t)(2,14)  = -0.030;  (*t)(2,14)  =   0.030;
1157       (*t)(0,15)  = 0.750;  (*t)(2,15)  = -0.030;  (*t)(2,15)  =   0.030;
1158       (*t)(0,16)  = 0.800;  (*t)(2,16)  = -0.030;  (*t)(2,16)  =   0.030;
1159       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.030;  (*t)(2,17)  =   0.030;
1160       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.030;  (*t)(2,18)  =   0.030;
1161       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.028;  (*t)(2,19)  =   0.028;
1162       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.028;  (*t)(2,20)  =   0.028;
1163       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.028;  (*t)(2,21)  =   0.028;
1164       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.026;  (*t)(2,22)  =   0.028;
1165       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.024;  (*t)(2,23)  =   0.028;
1166       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.020;  (*t)(2,24)  =   0.028;
1167       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.018;  (*t)(2,25)  =   0.028;
1168       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.016;  (*t)(2,26)  =   0.028;
1169       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.014;  (*t)(2,27)  =   0.028;
1170       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.012;  (*t)(2,28)  =   0.026;
1171       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.010;  (*t)(2,29)  =   0.026;
1172       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.008;  (*t)(2,30)  =   0.026;
1173       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.008;  (*t)(2,31)  =   0.024;
1174       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.006;  (*t)(2,32)  =   0.024;
1175       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.004;  (*t)(2,33)  =   0.024;
1176       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.004;  (*t)(2,34)  =   0.024;
1177       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.002;  (*t)(2,35)  =   0.024;
1178       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.002;  (*t)(2,36)  =   0.024;
1179       (*t)(0,37)  = 2.700;  (*t)(2,37)  = 0.000;  (*t)(2,37)  =   0.024;
1180       (*t)(0,38)  = 2.800;  (*t)(2,38)  = 0.000;  (*t)(2,38)  =   0.026;
1181       (*t)(0,39)  = 2.900;  (*t)(2,39)  = 0.000;  (*t)(2,39)  =   0.024;
1182       (*t)(0,40)  = 3.000;  (*t)(2,40)  = 0.002;  (*t)(2,40)  =   0.026;
1183       (*t)(0,41)  = 3.100;  (*t)(2,41)  = 0.002;  (*t)(2,41)  =   0.026;
1184       (*t)(0,42)  = 3.200;  (*t)(2,42)  = 0.002;  (*t)(2,42)  =   0.026;
1185       (*t)(0,43)  = 3.300;  (*t)(2,43)  = 0.002;  (*t)(2,43)  =   0.026;
1186       (*t)(0,44)  = 3.400;  (*t)(2,44)  = 0.002;  (*t)(2,44)  =   0.026;
1187       (*t)(0,45)  = 3.500;  (*t)(2,45)  = 0.002;  (*t)(2,45)  =   0.026;
1188       (*t)(0,46)  = 3.600;  (*t)(2,46)  = 0.002;  (*t)(2,46)  =   0.026;
1189       (*t)(0,47)  = 3.700;  (*t)(2,47)  = 0.002;  (*t)(2,47)  =   0.026;
1190       (*t)(0,48)  = 3.800;  (*t)(2,48)  = 0.002;  (*t)(2,48)  =   0.026;
1191       (*t)(0,49)  = 3.900;  (*t)(2,49)  = 0.004;  (*t)(2,49)  =   0.024;
1192       (*t)(0,50)  = 4.000;  (*t)(2,50)  = 0.004;  (*t)(2,50)  =   0.026;
1193       (*t)(0,51)  = 4.100;  (*t)(2,51)  = 0.004;  (*t)(2,51)  =   0.026;
1194       (*t)(0,52)  = 4.200;  (*t)(2,52)  = 0.004;  (*t)(2,52)  =   0.024;
1195       (*t)(0,53)  = 4.300;  (*t)(2,53)  = 0.006;  (*t)(2,53)  =   0.024;
1196       (*t)(0,54)  = 4.400;  (*t)(2,54)  = 0.000;  (*t)(2,54)  =   0.000;
1197       (*t)(0,55)  = 4.500;  (*t)(2,55)  = 0.000;  (*t)(2,55)  =   0.000;
1198       (*t)(0,56)  = 4.600;  (*t)(2,56)  = 0.000;  (*t)(2,56)  =   0.000;
1199       (*t)(0,57)  = 4.700;  (*t)(2,57)  = 0.000;  (*t)(2,57)  =   0.000;
1200       (*t)(0,58)  = 4.800;  (*t)(2,58)  = 0.000;  (*t)(2,58)  =   0.000;
1201       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1202       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000;
1203     }
1204     else
1205     if (fParticleID==AliPID::kProton)
1206     {
1207       //TOF protons, 0.9 purity
1208       t = new TMatrixF(3,61);
1209       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1210       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1211       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1212       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1213       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1214       (*t)(0,5)  = 0.250;  (*t)(2,5)  = 0.000;  (*t)(2,5)  =   0.000;
1215       (*t)(0,6)  = 0.300;  (*t)(2,6)  = 0.000;  (*t)(2,6)  =   0.000;
1216       (*t)(0,7)  = 0.350;  (*t)(2,7)  = 0.000;  (*t)(2,7)  =   0.000;
1217       (*t)(0,8)  = 0.400;  (*t)(2,8)  = 0.000;  (*t)(2,8)  =   0.000;
1218       (*t)(0,9)  = 0.450;  (*t)(2,9)  = 0.000;  (*t)(2,9)  =   0.000;
1219       (*t)(0,10)  = 0.500;  (*t)(2,10)  = 0.000;  (*t)(2,10)  =   0.000;
1220       (*t)(0,11)  = 0.550;  (*t)(2,11)  = 0.000;  (*t)(2,11)  =   0.000;
1221       (*t)(0,12)  = 0.600;  (*t)(2,12)  = 0.000;  (*t)(2,12)  =   0.000;
1222       (*t)(0,13)  = 0.650;  (*t)(2,13)  = 0.000;  (*t)(2,13)  =   0.000;
1223       (*t)(0,14)  = 0.700;  (*t)(2,14)  = 0.000;  (*t)(2,14)  =   0.000;
1224       (*t)(0,15)  = 0.750;  (*t)(2,15)  = 0.000;  (*t)(2,15)  =   0.000;
1225       (*t)(0,16)  = 0.800;  (*t)(2,16)  = 0.000;  (*t)(2,16)  =   0.000;
1226       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.070;  (*t)(2,17)  =   0.070;
1227       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.072;  (*t)(2,18)  =   0.072;
1228       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.072;  (*t)(2,19)  =   0.072;
1229       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.074;  (*t)(2,20)  =   0.074;
1230       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.032;  (*t)(2,21)  =   0.032;
1231       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.026;  (*t)(2,22)  =   0.026;
1232       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.026;  (*t)(2,23)  =   0.026;
1233       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.024;  (*t)(2,24)  =   0.024;
1234       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.024;  (*t)(2,25)  =   0.024;
1235       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.026;  (*t)(2,26)  =   0.026;
1236       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.026;  (*t)(2,27)  =   0.026;
1237       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.026;  (*t)(2,28)  =   0.026;
1238       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.026;  (*t)(2,29)  =   0.026;
1239       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.026;  (*t)(2,30)  =   0.026;
1240       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.026;  (*t)(2,31)  =   0.026;
1241       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.026;  (*t)(2,32)  =   0.024;
1242       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.028;  (*t)(2,33)  =   0.022;
1243       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.028;  (*t)(2,34)  =   0.020;
1244       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.028;  (*t)(2,35)  =   0.018;
1245       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.028;  (*t)(2,36)  =   0.016;
1246       (*t)(0,37)  = 2.700;  (*t)(2,37)  = -0.028;  (*t)(2,37)  =   0.016;
1247       (*t)(0,38)  = 2.800;  (*t)(2,38)  = -0.030;  (*t)(2,38)  =   0.014;
1248       (*t)(0,39)  = 2.900;  (*t)(2,39)  = -0.030;  (*t)(2,39)  =   0.012;
1249       (*t)(0,40)  = 3.000;  (*t)(2,40)  = -0.030;  (*t)(2,40)  =   0.012;
1250       (*t)(0,41)  = 3.100;  (*t)(2,41)  = -0.030;  (*t)(2,41)  =   0.010;
1251       (*t)(0,42)  = 3.200;  (*t)(2,42)  = -0.030;  (*t)(2,42)  =   0.010;
1252       (*t)(0,43)  = 3.300;  (*t)(2,43)  = -0.030;  (*t)(2,43)  =   0.010;
1253       (*t)(0,44)  = 3.400;  (*t)(2,44)  = -0.030;  (*t)(2,44)  =   0.008;
1254       (*t)(0,45)  = 3.500;  (*t)(2,45)  = -0.030;  (*t)(2,45)  =   0.008;
1255       (*t)(0,46)  = 3.600;  (*t)(2,46)  = -0.030;  (*t)(2,46)  =   0.008;
1256       (*t)(0,47)  = 3.700;  (*t)(2,47)  = -0.030;  (*t)(2,47)  =   0.006;
1257       (*t)(0,48)  = 3.800;  (*t)(2,48)  = -0.030;  (*t)(2,48)  =   0.006;
1258       (*t)(0,49)  = 3.900;  (*t)(2,49)  = -0.030;  (*t)(2,49)  =   0.006;
1259       (*t)(0,50)  = 4.000;  (*t)(2,50)  = -0.028;  (*t)(2,50)  =   0.004;
1260       (*t)(0,51)  = 4.100;  (*t)(2,51)  = -0.030;  (*t)(2,51)  =   0.004;
1261       (*t)(0,52)  = 4.200;  (*t)(2,52)  = -0.030;  (*t)(2,52)  =   0.004;
1262       (*t)(0,53)  = 4.300;  (*t)(2,53)  = -0.028;  (*t)(2,53)  =   0.002;
1263       (*t)(0,54)  = 4.400;  (*t)(2,54)  = -0.030;  (*t)(2,54)  =   0.002;
1264       (*t)(0,55)  = 4.500;  (*t)(2,55)  = -0.028;  (*t)(2,55)  =   0.002;
1265       (*t)(0,56)  = 4.600;  (*t)(2,56)  = -0.028;  (*t)(2,56)  =   0.002;
1266       (*t)(0,57)  = 4.700;  (*t)(2,57)  = -0.028;  (*t)(2,57)  =   0.000;
1267       (*t)(0,58)  = 4.800;  (*t)(2,58)  = -0.028;  (*t)(2,58)  =   0.002;
1268       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1269       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000;
1270     }
1271     else
1272     if (fParticleID==AliPID::kKaon)
1273     {
1274       //TOF kaons, 0.9 purity
1275       t = new TMatrixF(3,61);
1276       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1277       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1278       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1279       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1280       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1281       (*t)(0,5)  = 0.250;  (*t)(2,5)  = 0.000;  (*t)(2,5)  =   0.000;
1282       (*t)(0,6)  = 0.300;  (*t)(2,6)  = 0.000;  (*t)(2,6)  =   0.000;
1283       (*t)(0,7)  = 0.350;  (*t)(2,7)  = 0.000;  (*t)(2,7)  =   0.000;
1284       (*t)(0,8)  = 0.400;  (*t)(2,8)  = 0.000;  (*t)(2,8)  =   0.000;
1285       (*t)(0,9)  = 0.450;  (*t)(2,9)  = 0.000;  (*t)(2,9)  =   0.000;
1286       (*t)(0,10)  = 0.500;  (*t)(2,10)  = 0.000;  (*t)(2,10)  =   0.000;
1287       (*t)(0,11)  = 0.550;  (*t)(2,11)  = -0.026;  (*t)(2,11)  =   0.026;
1288       (*t)(0,12)  = 0.600;  (*t)(2,12)  = -0.026;  (*t)(2,12)  =   0.026;
1289       (*t)(0,13)  = 0.650;  (*t)(2,13)  = -0.026;  (*t)(2,13)  =   0.026;
1290       (*t)(0,14)  = 0.700;  (*t)(2,14)  = -0.026;  (*t)(2,14)  =   0.026;
1291       (*t)(0,15)  = 0.750;  (*t)(2,15)  = -0.026;  (*t)(2,15)  =   0.026;
1292       (*t)(0,16)  = 0.800;  (*t)(2,16)  = -0.026;  (*t)(2,16)  =   0.026;
1293       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.024;  (*t)(2,17)  =   0.024;
1294       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.024;  (*t)(2,18)  =   0.024;
1295       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.024;  (*t)(2,19)  =   0.024;
1296       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.024;  (*t)(2,20)  =   0.024;
1297       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.024;  (*t)(2,21)  =   0.024;
1298       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.024;  (*t)(2,22)  =   0.022;
1299       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.024;  (*t)(2,23)  =   0.020;
1300       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.026;  (*t)(2,24)  =   0.016;
1301       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.028;  (*t)(2,25)  =   0.014;
1302       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.028;  (*t)(2,26)  =   0.012;
1303       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.028;  (*t)(2,27)  =   0.010;
1304       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.028;  (*t)(2,28)  =   0.010;
1305       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.028;  (*t)(2,29)  =   0.008;
1306       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.028;  (*t)(2,30)  =   0.006;
1307       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.026;  (*t)(2,31)  =   0.006;
1308       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.024;  (*t)(2,32)  =   0.004;
1309       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.020;  (*t)(2,33)  =   0.002;
1310       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.020;  (*t)(2,34)  =   0.002;
1311       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.018;  (*t)(2,35)  =   0.000;
1312       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.016;  (*t)(2,36)  =   0.000;
1313       (*t)(0,37)  = 2.700;  (*t)(2,37)  = -0.014;  (*t)(2,37)  =   -0.002;
1314       (*t)(0,38)  = 2.800;  (*t)(2,38)  = -0.014;  (*t)(2,38)  =   -0.004;
1315       (*t)(0,39)  = 2.900;  (*t)(2,39)  = -0.012;  (*t)(2,39)  =   -0.004;
1316       (*t)(0,40)  = 3.000;  (*t)(2,40)  = -0.010;  (*t)(2,40)  =   -0.006;
1317       (*t)(0,41)  = 3.100;  (*t)(2,41)  = 0.000;  (*t)(2,41)  =   0.000;
1318       (*t)(0,42)  = 3.200;  (*t)(2,42)  = 0.000;  (*t)(2,42)  =   0.000;
1319       (*t)(0,43)  = 3.300;  (*t)(2,43)  = 0.000;  (*t)(2,43)  =   0.000;
1320       (*t)(0,44)  = 3.400;  (*t)(2,44)  = 0.000;  (*t)(2,44)  =   0.000;
1321       (*t)(0,45)  = 3.500;  (*t)(2,45)  = 0.000;  (*t)(2,45)  =   0.000;
1322       (*t)(0,46)  = 3.600;  (*t)(2,46)  = 0.000;  (*t)(2,46)  =   0.000;
1323       (*t)(0,47)  = 3.700;  (*t)(2,47)  = 0.000;  (*t)(2,47)  =   0.000;
1324       (*t)(0,48)  = 3.800;  (*t)(2,48)  = 0.000;  (*t)(2,48)  =   0.000;
1325       (*t)(0,49)  = 3.900;  (*t)(2,49)  = 0.000;  (*t)(2,49)  =   0.000;
1326       (*t)(0,50)  = 4.000;  (*t)(2,50)  = 0.000;  (*t)(2,50)  =   0.000;
1327       (*t)(0,51)  = 4.100;  (*t)(2,51)  = 0.000;  (*t)(2,51)  =   0.000;
1328       (*t)(0,52)  = 4.200;  (*t)(2,52)  = 0.000;  (*t)(2,52)  =   0.000;
1329       (*t)(0,53)  = 4.300;  (*t)(2,53)  = 0.000;  (*t)(2,53)  =   0.000;
1330       (*t)(0,54)  = 4.400;  (*t)(2,54)  = 0.000;  (*t)(2,54)  =   0.000;
1331       (*t)(0,55)  = 4.500;  (*t)(2,55)  = 0.000;  (*t)(2,55)  =   0.000;
1332       (*t)(0,56)  = 4.600;  (*t)(2,56)  = 0.000;  (*t)(2,56)  =   0.000;
1333       (*t)(0,57)  = 4.700;  (*t)(2,57)  = 0.000;  (*t)(2,57)  =   0.000;
1334       (*t)(0,58)  = 4.800;  (*t)(2,58)  = 0.000;  (*t)(2,58)  =   0.000;
1335       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1336       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000;
1337     }
1338     fTOFpidCuts=t;
1339   }
1340 }
1341
1342 //-----------------------------------------------------------------------
1343 // part added by F. Noferini (some methods)
1344 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
1345
1346   Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1347
1348   if (! goodtrack)
1349        return kFALSE;
1350
1351   Int_t pdg = GetESDPdg(track,"bayesianALL");
1352   //  printf("pdg set to %i\n",pdg);
1353
1354   Int_t pid = 0;
1355   Float_t prob = 0;
1356   switch (fParticleID)
1357   {
1358     case AliPID::kPion:
1359       pid=211;
1360       prob = fProbBayes[2];
1361       break;
1362     case AliPID::kKaon:
1363       pid=321;
1364       prob = fProbBayes[3];
1365      break;
1366     case AliPID::kProton:
1367       pid=2212;
1368       prob = fProbBayes[4];
1369       break;
1370     case AliPID::kElectron:
1371       pid=-11;
1372        prob = fProbBayes[0];
1373      break;
1374     default:
1375       return kFALSE;
1376   }
1377
1378   //  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);
1379   if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1380     if(!fCutCharge)
1381       return kTRUE;
1382     else if (fCutCharge && fCharge * track->GetSign() > 0)
1383       return kTRUE;
1384   }
1385   return kFALSE;
1386 }
1387 //-----------------------------------------------------------------------
1388 Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1389   Int_t pdg = 0;
1390   Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1391   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1392
1393   if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1394     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1395     Double_t rcc=0.;
1396
1397     Float_t pt = track->Pt();
1398
1399     Int_t iptesd = 0;
1400     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1401
1402     if(cPi < 0){
1403       c[0] = fC[iptesd][0];
1404       c[1] = fC[iptesd][1];
1405       c[2] = fC[iptesd][2];
1406       c[3] = fC[iptesd][3];
1407       c[4] = fC[iptesd][4];
1408     }
1409     else{
1410       c[0] = 0.0;
1411       c[1] = 0.0;
1412       c[2] = cPi;
1413       c[3] = cKa;
1414       c[4] = cPr;
1415     }
1416
1417     Double_t r1[10]; track->GetTOFpid(r1);
1418
1419     Int_t i;
1420     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1421
1422     Double_t w[10];
1423     for (i=0; i<5; i++){
1424         w[i]=c[i]*r1[i]/rcc;
1425         fProbBayes[i] = w[i];
1426     }
1427     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1428       pdg = 211*Int_t(track->GetSign());
1429     }
1430     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1431       pdg = 2212*Int_t(track->GetSign());
1432     }
1433     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1434       pdg = 321*Int_t(track->GetSign());
1435     }
1436     else if (w[0]>=w[1]) { //electrons
1437       pdg = -11*Int_t(track->GetSign());
1438     }
1439     else{ // muon
1440       pdg = -13*Int_t(track->GetSign());
1441     }
1442   }
1443
1444   else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1445     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1446     Double_t rcc=0.;
1447
1448     Float_t pt = track->Pt();
1449
1450     Int_t iptesd = 0;
1451     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1452
1453     if(cPi < 0){
1454       c[0] = fC[iptesd][0];
1455       c[1] = fC[iptesd][1];
1456       c[2] = fC[iptesd][2];
1457       c[3] = fC[iptesd][3];
1458       c[4] = fC[iptesd][4];
1459     }
1460     else{
1461       c[0] = 0.0;
1462       c[1] = 0.0;
1463       c[2] = cPi;
1464       c[3] = cKa;
1465       c[4] = cPr;
1466     }
1467
1468     Double_t r1[10]; track->GetTPCpid(r1);
1469
1470     Int_t i;
1471     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1472
1473     Double_t w[10];
1474     for (i=0; i<5; i++){
1475         w[i]=c[i]*r1[i]/rcc;
1476         fProbBayes[i] = w[i];
1477     }
1478     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1479       pdg = 211*Int_t(track->GetSign());
1480     }
1481     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1482       pdg = 2212*Int_t(track->GetSign());
1483     }
1484     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1485       pdg = 321*Int_t(track->GetSign());
1486     }
1487     else if (w[0]>=w[1]) { //electrons
1488       pdg = -11*Int_t(track->GetSign());
1489     }
1490     else{ // muon
1491       pdg = -13*Int_t(track->GetSign());
1492     }
1493   }
1494
1495   else if(strstr(option,"bayesianALL")){
1496     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1497     Double_t rcc=0.;
1498
1499     Float_t pt = track->Pt();
1500
1501     Int_t iptesd = 0;
1502     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1503
1504     if(cPi < 0){
1505       c[0] = fC[iptesd][0];
1506       c[1] = fC[iptesd][1];
1507       c[2] = fC[iptesd][2];
1508       c[3] = fC[iptesd][3];
1509       c[4] = fC[iptesd][4];
1510     }
1511     else{
1512       c[0] = 0.0;
1513       c[1] = 0.0;
1514       c[2] = cPi;
1515       c[3] = cKa;
1516       c[4] = cPr;
1517     }
1518
1519     Double_t r1[10]; track->GetTOFpid(r1);
1520     Double_t r2[10]; track->GetTPCpid(r2);
1521
1522     Int_t i;
1523     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1524
1525
1526     Double_t w[10];
1527     for (i=0; i<5; i++){
1528         w[i]=c[i]*r1[i]*r2[i]/rcc;
1529         fProbBayes[i] = w[i];
1530     }
1531
1532     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1533       pdg = 211*Int_t(track->GetSign());
1534     }
1535     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1536       pdg = 2212*Int_t(track->GetSign());
1537     }
1538     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1539       pdg = 321*Int_t(track->GetSign());
1540     }
1541     else if (w[0]>=w[1]) { //electrons
1542       pdg = -11*Int_t(track->GetSign());
1543     }
1544     else{ // muon
1545       pdg = -13*Int_t(track->GetSign());
1546     }
1547   }
1548
1549   else if(strstr(option,"sigmacutTOF")){
1550     printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1551     Float_t p = track->P();
1552
1553     // Take expected times
1554     Double_t exptimes[5];
1555     track->GetIntegratedTimes(exptimes);
1556
1557     // Take resolution for TOF response
1558     // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1559     Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1560
1561     if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1562       pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1563     }
1564   }
1565
1566   else{
1567     printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1568   }
1569
1570   return pdg;
1571 }
1572 //-----------------------------------------------------------------------
1573 void AliFlowTrackCuts::SetPriors(){
1574   // set abbundancies
1575     fBinLimitPID[0] = 0.30;
1576     fC[0][0] = 0.015;
1577     fC[0][1] = 0.015;
1578     fC[0][2] = 1;
1579     fC[0][3] = 0.0025;
1580     fC[0][4] = 0.000015;
1581     fBinLimitPID[1] = 0.35;
1582     fC[1][0] = 0.015;
1583     fC[1][1] = 0.015;
1584     fC[1][2] = 1;
1585     fC[1][3] = 0.01;
1586     fC[1][4] = 0.001;
1587     fBinLimitPID[2] = 0.40;
1588     fC[2][0] = 0.015;
1589     fC[2][1] = 0.015;
1590     fC[2][2] = 1;
1591     fC[2][3] = 0.026;
1592     fC[2][4] = 0.004;
1593     fBinLimitPID[3] = 0.45;
1594     fC[3][0] = 0.015;
1595     fC[3][1] = 0.015;
1596     fC[3][2] = 1;
1597     fC[3][3] = 0.026;
1598     fC[3][4] = 0.004;
1599     fBinLimitPID[4] = 0.50;
1600     fC[4][0] = 0.015;
1601     fC[4][1] = 0.015;
1602     fC[4][2] = 1.000000;
1603     fC[4][3] = 0.05;
1604     fC[4][4] = 0.01;
1605     fBinLimitPID[5] = 0.60;
1606     fC[5][0] = 0.012;
1607     fC[5][1] = 0.012;
1608     fC[5][2] = 1;
1609     fC[5][3] = 0.085;
1610     fC[5][4] = 0.022;
1611     fBinLimitPID[6] = 0.70;
1612     fC[6][0] = 0.01;
1613     fC[6][1] = 0.01;
1614     fC[6][2] = 1;
1615     fC[6][3] = 0.12;
1616     fC[6][4] = 0.036;
1617     fBinLimitPID[7] = 0.80;
1618     fC[7][0] = 0.0095;
1619     fC[7][1] = 0.0095;
1620     fC[7][2] = 1;
1621     fC[7][3] = 0.15;
1622     fC[7][4] = 0.05;
1623     fBinLimitPID[8] = 0.90;
1624     fC[8][0] = 0.0085;
1625     fC[8][1] = 0.0085;
1626     fC[8][2] = 1;
1627     fC[8][3] = 0.18;
1628     fC[8][4] = 0.074;
1629     fBinLimitPID[9] = 1;
1630     fC[9][0] = 0.008;
1631     fC[9][1] = 0.008;
1632     fC[9][2] = 1;
1633     fC[9][3] = 0.22;
1634     fC[9][4] = 0.1;
1635     fBinLimitPID[10] = 1.20;
1636     fC[10][0] = 0.007;
1637     fC[10][1] = 0.007;
1638     fC[10][2] = 1;
1639     fC[10][3] = 0.28;
1640     fC[10][4] = 0.16;
1641     fBinLimitPID[11] = 1.40;
1642     fC[11][0] = 0.0066;
1643     fC[11][1] = 0.0066;
1644     fC[11][2] = 1;
1645     fC[11][3] = 0.35;
1646     fC[11][4] = 0.23;
1647     fBinLimitPID[12] = 1.60;
1648     fC[12][0] = 0.0075;
1649     fC[12][1] = 0.0075;
1650     fC[12][2] = 1;
1651     fC[12][3] = 0.40;
1652     fC[12][4] = 0.31;
1653     fBinLimitPID[13] = 1.80;
1654     fC[13][0] = 0.0062;
1655     fC[13][1] = 0.0062;
1656     fC[13][2] = 1;
1657     fC[13][3] = 0.45;
1658     fC[13][4] = 0.39;
1659     fBinLimitPID[14] = 2.00;
1660     fC[14][0] = 0.005;
1661     fC[14][1] = 0.005;
1662     fC[14][2] = 1;
1663     fC[14][3] = 0.46;
1664     fC[14][4] = 0.47;
1665     fBinLimitPID[15] = 2.20;
1666     fC[15][0] = 0.0042;
1667     fC[15][1] = 0.0042;
1668     fC[15][2] = 1;
1669     fC[15][3] = 0.5;
1670     fC[15][4] = 0.55;
1671     fBinLimitPID[16] = 2.40;
1672     fC[16][0] = 0.007;
1673     fC[16][1] = 0.007;
1674     fC[16][2] = 1;
1675     fC[16][3] = 0.5;
1676     fC[16][4] = 0.6;
1677
1678     for(Int_t i=17;i<fnPIDptBin;i++){
1679         fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1680         fC[i][0] = fC[13][0];
1681         fC[i][1] = fC[13][1];
1682         fC[i][2] = fC[13][2];
1683         fC[i][3] = fC[13][3];
1684         fC[i][4] = fC[13][4];
1685     }
1686 }
1687 // end part added by F. Noferini
1688 //-----------------------------------------------------------------------
1689
1690
1691 //-----------------------------------------------------------------------
1692 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
1693 {
1694   //get the name of the particle id source
1695   switch (s)
1696   {
1697     case kTPCdedx:
1698       return "TPCdedx";
1699     case kTOFbeta:
1700       return "TOFbeta";
1701     case kTPCpid:
1702       return "TPCpid";
1703     case kTOFpid:
1704       return "TOFpid";
1705     case kTOFbayesian:
1706       return "TOFbayesianPID";
1707     default:
1708       return "NOPID";
1709   }
1710 }
1711
1712 //-----------------------------------------------------------------------
1713 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
1714 {
1715   //return the name of the selected parameter type
1716   switch (type)
1717   {
1718     case kMC:
1719       return "MC";
1720     case kGlobal:
1721       return "ESD global";
1722     case kESD_TPConly:
1723       return "TPC only";
1724     case kESD_SPDtracklet:
1725       return "SPD tracklet";
1726     default:
1727       return "unknown";
1728   }
1729 }
1730