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