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