PWG/FLOW/Tasks/AliFlowTrackCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **********************************************************i****************/
  15
  16 /* $Id$ */
  17
  18 // AliFlowTrackCuts:
  19 // Data selection for flow framework
  20 //
  21 // origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
  22 // mods:   Redmer A. Bertens (rbertens@cern.ch)
  23 //
  24 // This class gurantees consistency of cut methods, trackparameter
  25 // selection (global tracks, TPC only, etc..) and parameter mixing
  26 // in the flow framework. Transparently handles different input types:
  27 // ESD, MC, AOD.
  28 // This class works in 2 steps: first the requested track parameters are
  29 // constructed (to be set by SetParamType() ), then cuts are applied.
  30 // the constructed track can be requested AFTER checking the cuts by
  31 // calling GetTrack(), in this case the cut object stays in control,
  32 // caller does not have to delete the track.
  33 // Additionally caller can request an AliFlowTrack object to be constructed
  34 // according the parameter mixing scenario requested by SetParamMix().
  35 // AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
  36 // so caller needs to take care of the freshly created object.
  37
  38 #include <limits.h>
  39 #include <float.h>
  40 #include <TMatrix.h>
  41 #include "TMCProcess.h"
  42 #include "TParticle.h"
  43 #include "TH2F.h"
  44 #include "AliStack.h"
  45 #include "TBrowser.h"
  46 #include "AliMCEvent.h"
  47 #include "AliESDEvent.h"
  48 #include "AliAODEvent.h"
  49 #include "AliVParticle.h"
  50 #include "AliVVZERO.h"
  51 #include "AliMCParticle.h"
  52 #include "AliESDkink.h"
  53 #include "AliESDv0.h"
  54 #include "AliESDtrack.h"
  55 #include "AliESDMuonTrack.h"   // XZhang 20120604
  56 #include "AliMultiplicity.h"
  57 #include "AliAODTrack.h"
  58 #include "AliAODTracklets.h"   // XZhang 20120615
  59 #include "AliFlowTrackSimple.h"
  60 #include "AliFlowTrack.h"
  61 #include "AliFlowTrackCuts.h"
  62 #include "AliLog.h"
  63 #include "AliESDpid.h"
  64 #include "AliESDPmdTrack.h"
  65 #include "AliESDUtils.h" //TODO
  66 #include "AliFlowBayesianPID.h"
  67 #include "AliFlowCandidateTrack.h"
  68 #include "AliKFParticle.h"
  69 #include "AliESDVZERO.h"
  70 #include "AliFlowCommonConstants.h"
  71
  72 ClassImp(AliFlowTrackCuts)
  73
  74 //-----------------------------------------------------------------------
  75 AliFlowTrackCuts::AliFlowTrackCuts():
  76   AliFlowTrackSimpleCuts(),
  77   fAliESDtrackCuts(NULL),
  78   fMuonTrackCuts(NULL),  // XZhang 20120604
  79   fQA(NULL),
  80   fCutMC(kFALSE),
  81   fCutMChasTrackReferences(kFALSE),
  82   fCutMCprocessType(kFALSE),
  83   fMCprocessType(kPNoProcess),
  84   fCutMCPID(kFALSE),
  85   fMCPID(0),
  86   fCutMCfirstMotherPID(kFALSE),
  87   fMCfirstMotherPID(0),
  88   fIgnoreSignInMCPID(kFALSE),
  89   fCutMCisPrimary(kFALSE),
  90   fRequireTransportBitForPrimaries(kTRUE),
  91   fMCisPrimary(kFALSE),
  92   fRequireCharge(kFALSE),
  93   fFakesAreOK(kTRUE),
  94   fCutSPDtrackletDeltaPhi(kFALSE),
  95   fSPDtrackletDeltaPhiMax(FLT_MAX),
  96   fSPDtrackletDeltaPhiMin(-FLT_MAX),
  97   fIgnoreTPCzRange(kFALSE),
  98   fIgnoreTPCzRangeMax(FLT_MAX),
  99   fIgnoreTPCzRangeMin(-FLT_MAX),
 100   fCutChi2PerClusterTPC(kFALSE),
 101   fMaxChi2PerClusterTPC(FLT_MAX),
 102   fMinChi2PerClusterTPC(-FLT_MAX),
 103   fCutNClustersTPC(kFALSE),
 104   fNClustersTPCMax(INT_MAX),
 105   fNClustersTPCMin(INT_MIN),
 106   fCutNClustersITS(kFALSE),
 107   fNClustersITSMax(INT_MAX),
 108   fNClustersITSMin(INT_MIN),
 109   fUseAODFilterBit(kTRUE),
 110   fAODFilterBit(1),
 111   fCutDCAToVertexXY(kFALSE),
 112   fCutDCAToVertexZ(kFALSE),
 113   fCutMinimalTPCdedx(kFALSE),
 114   fMinimalTPCdedx(0.),
 115   fLinearizeVZEROresponse(kFALSE),
 116   fCutPmdDet(kFALSE),
 117   fPmdDet(0),
 118   fCutPmdAdc(kFALSE),
 119   fPmdAdc(0.),
 120   fCutPmdNcell(kFALSE),
 121   fPmdNcell(0.),
 122   fMinKinkAngle(TMath::DegToRad()*2.),
 123   fMinKinkRadius(130.),
 124   fMaxKinkRadius(200.),
 125   fMinKinkQt(.05),
 126   fMaxKinkQt(.5),
 127   fMinKinkInvMassKmu(0.),
 128   fMaxKinkInvMassKmu(0.6),
 129   fForceTPCstandalone(kFALSE),
 130   fRequireKinkDaughters(kFALSE),
 131   fParamType(kGlobal),
 132   fParamMix(kPure),
 133   fKink(NULL),
 134   fV0(NULL),
 135   fTrack(NULL),
 136   fTrackMass(0.),
 137   fTrackPt(0.),
 138   fTrackPhi(0.),
 139   fTrackEta(0.),
 140   fTrackWeight(1.),
 141   fTrackLabel(INT_MIN),
 142   fMCevent(NULL),
 143   fMCparticle(NULL),
 144   fEvent(NULL),
 145   fTPCtrack(),
 146   fESDpid(),
 147   fBayesianResponse(NULL),
 148   fPIDsource(kTOFpid),
 149   fTPCpidCuts(NULL),
 150   fTOFpidCuts(NULL),
 151   fParticleID(AliPID::kUnknown),
 152   fParticleProbability(.9),
 153   fAllowTOFmismatchFlag(kFALSE),
 154   fRequireStrictTOFTPCagreement(kFALSE),
 155   fCutRejectElectronsWithTPCpid(kFALSE),
 156   fProbBayes(0.0),
 157   fCurrCentr(0.0),
 158   fVZEROgainEqualization(NULL),
 159   fApplyRecentering(kFALSE),
 160   fVZEROgainEqualizationPerRing(kFALSE)
 161 {
 162   //io constructor
 163   SetPriors(); //init arrays
 164
 165   // New PID procedure (Bayesian Combined PID)
 166   // allocating here is necessary because we don't
 167   // stream this member
 168   // TODO: fix streaming problems AliFlowBayesianPID
 169   fBayesianResponse = new AliFlowBayesianPID();
 170   fBayesianResponse->SetNewTrackParam();
 171   for(Int_t i(0); i < 4; i++) {
 172       fVZEROApol[i] = 0;
 173       fVZEROCpol[i] = 0;
 174   }
 175   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
 176
 177 }
 178
 179 //-----------------------------------------------------------------------
 180 AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
 181   AliFlowTrackSimpleCuts(name),
 182   fAliESDtrackCuts(NULL),
 183   fMuonTrackCuts(NULL),  // XZhang 20120604
 184   fQA(NULL),
 185   fCutMC(kFALSE),
 186   fCutMChasTrackReferences(kFALSE),
 187   fCutMCprocessType(kFALSE),
 188   fMCprocessType(kPNoProcess),
 189   fCutMCPID(kFALSE),
 190   fMCPID(0),
 191   fCutMCfirstMotherPID(kFALSE),
 192   fMCfirstMotherPID(0),
 193   fIgnoreSignInMCPID(kFALSE),
 194   fCutMCisPrimary(kFALSE),
 195   fRequireTransportBitForPrimaries(kTRUE),
 196   fMCisPrimary(kFALSE),
 197   fRequireCharge(kFALSE),
 198   fFakesAreOK(kTRUE),
 199   fCutSPDtrackletDeltaPhi(kFALSE),
 200   fSPDtrackletDeltaPhiMax(FLT_MAX),
 201   fSPDtrackletDeltaPhiMin(-FLT_MAX),
 202   fIgnoreTPCzRange(kFALSE),
 203   fIgnoreTPCzRangeMax(FLT_MAX),
 204   fIgnoreTPCzRangeMin(-FLT_MAX),
 205   fCutChi2PerClusterTPC(kFALSE),
 206   fMaxChi2PerClusterTPC(FLT_MAX),
 207   fMinChi2PerClusterTPC(-FLT_MAX),
 208   fCutNClustersTPC(kFALSE),
 209   fNClustersTPCMax(INT_MAX),
 210   fNClustersTPCMin(INT_MIN),
 211   fCutNClustersITS(kFALSE),
 212   fNClustersITSMax(INT_MAX),
 213   fNClustersITSMin(INT_MIN),
 214   fUseAODFilterBit(kTRUE),
 215   fAODFilterBit(1),
 216   fCutDCAToVertexXY(kFALSE),
 217   fCutDCAToVertexZ(kFALSE),
 218   fCutMinimalTPCdedx(kFALSE),
 219   fMinimalTPCdedx(0.),
 220   fLinearizeVZEROresponse(kFALSE),
 221   fCutPmdDet(kFALSE),
 222   fPmdDet(0),
 223   fCutPmdAdc(kFALSE),
 224   fPmdAdc(0.),
 225   fCutPmdNcell(kFALSE),
 226   fPmdNcell(0.),
 227   fMinKinkAngle(TMath::DegToRad()*2.),
 228   fMinKinkRadius(130.),
 229   fMaxKinkRadius(200.),
 230   fMinKinkQt(0.05),
 231   fMaxKinkQt(0.5),
 232   fMinKinkInvMassKmu(0.0),
 233   fMaxKinkInvMassKmu(0.6),
 234   fForceTPCstandalone(kFALSE),
 235   fRequireKinkDaughters(kFALSE),
 236   fParamType(kGlobal),
 237   fParamMix(kPure),
 238   fKink(NULL),
 239   fV0(NULL),
 240   fTrack(NULL),
 241   fTrackMass(0.),
 242   fTrackPt(0.),
 243   fTrackPhi(0.),
 244   fTrackEta(0.),
 245   fTrackWeight(1.),
 246   fTrackLabel(INT_MIN),
 247   fMCevent(NULL),
 248   fMCparticle(NULL),
 249   fEvent(NULL),
 250   fTPCtrack(),
 251   fESDpid(),
 252   fBayesianResponse(NULL),
 253   fPIDsource(kTOFpid),
 254   fTPCpidCuts(NULL),
 255   fTOFpidCuts(NULL),
 256   fParticleID(AliPID::kUnknown),
 257   fParticleProbability(.9),
 258   fAllowTOFmismatchFlag(kFALSE),
 259   fRequireStrictTOFTPCagreement(kFALSE),
 260   fCutRejectElectronsWithTPCpid(kFALSE),
 261   fProbBayes(0.0),
 262   fCurrCentr(0.0),
 263   fVZEROgainEqualization(NULL),
 264   fApplyRecentering(kFALSE),
 265   fVZEROgainEqualizationPerRing(kFALSE)
 266 {
 267   //constructor
 268   SetTitle("AliFlowTrackCuts");
 269   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
 270                                                     2.63394e+01,
 271                                                     5.04114e-11,
 272                                                     2.12543e+00,
 273                                                     4.88663e+00 );
 274   SetPriors(); //init arrays
 275
 276   // New PID procedure (Bayesian Combined PID)
 277   fBayesianResponse = new AliFlowBayesianPID();
 278   fBayesianResponse->SetNewTrackParam();
 279   for(Int_t i(0); i < 4; i++) {
 280       fVZEROApol[i] = 0;
 281       fVZEROCpol[i] = 0;
 282   }
 283   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
 284 }
 285
 286 //-----------------------------------------------------------------------
 287 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
 288   AliFlowTrackSimpleCuts(that),
 289   fAliESDtrackCuts(NULL),
 290   fMuonTrackCuts(NULL),  // XZhang 20120604
 291   fQA(NULL),
 292   fCutMC(that.fCutMC),
 293   fCutMChasTrackReferences(that.fCutMChasTrackReferences),
 294   fCutMCprocessType(that.fCutMCprocessType),
 295   fMCprocessType(that.fMCprocessType),
 296   fCutMCPID(that.fCutMCPID),
 297   fMCPID(that.fMCPID),
 298   fCutMCfirstMotherPID(that.fCutMCfirstMotherPID),
 299   fMCfirstMotherPID(that.fMCfirstMotherPID),
 300   fIgnoreSignInMCPID(that.fIgnoreSignInMCPID),
 301   fCutMCisPrimary(that.fCutMCisPrimary),
 302   fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
 303   fMCisPrimary(that.fMCisPrimary),
 304   fRequireCharge(that.fRequireCharge),
 305   fFakesAreOK(that.fFakesAreOK),
 306   fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
 307   fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
 308   fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
 309   fIgnoreTPCzRange(that.fIgnoreTPCzRange),
 310   fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
 311   fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
 312   fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
 313   fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
 314   fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
 315   fCutNClustersTPC(that.fCutNClustersTPC),
 316   fNClustersTPCMax(that.fNClustersTPCMax),
 317   fNClustersTPCMin(that.fNClustersTPCMin),
 318   fCutNClustersITS(that.fCutNClustersITS),
 319   fNClustersITSMax(that.fNClustersITSMax),
 320   fNClustersITSMin(that.fNClustersITSMin),
 321   fUseAODFilterBit(that.fUseAODFilterBit),
 322   fAODFilterBit(that.fAODFilterBit),
 323   fCutDCAToVertexXY(that.fCutDCAToVertexXY),
 324   fCutDCAToVertexZ(that.fCutDCAToVertexZ),
 325   fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
 326   fMinimalTPCdedx(that.fMinimalTPCdedx),
 327   fLinearizeVZEROresponse(that.fLinearizeVZEROresponse),
 328   fCutPmdDet(that.fCutPmdDet),
 329   fPmdDet(that.fPmdDet),
 330   fCutPmdAdc(that.fCutPmdAdc),
 331   fPmdAdc(that.fPmdAdc),
 332   fCutPmdNcell(that.fCutPmdNcell),
 333   fPmdNcell(that.fPmdNcell),
 334   fMinKinkAngle(that.fMinKinkAngle),
 335   fMinKinkRadius(that.fMinKinkRadius),
 336   fMaxKinkRadius(that.fMaxKinkRadius),
 337   fMinKinkQt(that.fMinKinkQt),
 338   fMaxKinkQt(that.fMaxKinkQt),
 339   fMinKinkInvMassKmu(that.fMinKinkInvMassKmu),
 340   fMaxKinkInvMassKmu(that.fMaxKinkInvMassKmu),
 341   fForceTPCstandalone(that.fForceTPCstandalone),
 342   fRequireKinkDaughters(that.fRequireKinkDaughters),
 343   fParamType(that.fParamType),
 344   fParamMix(that.fParamMix),
 345   fKink(NULL),
 346   fV0(NULL),
 347   fTrack(NULL),
 348   fTrackMass(0.),
 349   fTrackPt(0.),
 350   fTrackPhi(0.),
 351   fTrackEta(0.),
 352   fTrackWeight(1.),
 353   fTrackLabel(INT_MIN),
 354   fMCevent(NULL),
 355   fMCparticle(NULL),
 356   fEvent(NULL),
 357   fTPCtrack(),
 358   fESDpid(that.fESDpid),
 359   fBayesianResponse(NULL),
 360   fPIDsource(that.fPIDsource),
 361   fTPCpidCuts(NULL),
 362   fTOFpidCuts(NULL),
 363   fParticleID(that.fParticleID),
 364   fParticleProbability(that.fParticleProbability),
 365   fAllowTOFmismatchFlag(that.fAllowTOFmismatchFlag),
 366   fRequireStrictTOFTPCagreement(that.fRequireStrictTOFTPCagreement),
 367   fCutRejectElectronsWithTPCpid(that.fCutRejectElectronsWithTPCpid),
 368   fProbBayes(0.0),
 369   fCurrCentr(0.0),
 370   fVZEROgainEqualization(NULL),
 371   fApplyRecentering(that.fApplyRecentering),
 372   fVZEROgainEqualizationPerRing(that.fVZEROgainEqualizationPerRing)
 373 {
 374   //copy constructor
 375   printf(" \n\n claling copy ctor \n\n" );
 376   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 377   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 378   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
 379   if (that.fMuonTrackCuts)   fMuonTrackCuts   = new AliMuonTrackCuts(*(that.fMuonTrackCuts));  // XZhang 20120604
 380   SetPriors(); //init arrays
 381   if (that.fQA) DefineHistograms();
 382
 383   // New PID procedure (Bayesian Combined PID)
 384   fBayesianResponse = new AliFlowBayesianPID();
 385   fBayesianResponse->SetNewTrackParam();
 386
 387   // VZERO gain calibration
 388   // no reason to init fVZEROgainEqualizationPerRing, will be initialized on node if necessary
 389   // pointer is set to NULL in initialization list of this constructor
 390 //  if (that.fVZEROgainEqualization) fVZEROgainEqualization = new TH1(*(that.fVZEROgainEqualization));
 391   for(Int_t i(0); i < 4; i++) { // no use to copy these guys since they're only initialized on worked node
 392       fVZEROApol[i] = 0.;
 393       fVZEROCpol[i] = 0.;
 394   }
 395   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = that.fUseVZERORing[i];
 396
 397 }
 398
 399 //-----------------------------------------------------------------------
 400 AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
 401 {
 402   //assignment
 403   if (this==&that) return *this;
 404
 405   AliFlowTrackSimpleCuts::operator=(that);
 406   //the following may seem excessive but if AliESDtrackCuts properly does copy and clone
 407   //this approach is better memory-fragmentation-wise in some cases
 408   if (that.fAliESDtrackCuts && fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
 409   if (that.fAliESDtrackCuts && !fAliESDtrackCuts) fAliESDtrackCuts=new AliESDtrackCuts(*(that.fAliESDtrackCuts));
 410   if (!that.fAliESDtrackCuts) delete fAliESDtrackCuts; fAliESDtrackCuts=NULL;
 411
 412   if ( that.fMuonTrackCuts &&  fMuonTrackCuts) *fMuonTrackCuts = *(that.fMuonTrackCuts);                   // XZhang 20120604
 413   if ( that.fMuonTrackCuts && !fMuonTrackCuts)  fMuonTrackCuts = new AliMuonTrackCuts(*(that.fMuonTrackCuts));  // XZhang 20120604
 414   if (!that.fMuonTrackCuts) delete fMuonTrackCuts; fMuonTrackCuts = NULL;                                  // XZhang 20120604
 415   if (!that.fVZEROgainEqualization) delete fVZEROgainEqualization; fVZEROgainEqualization = NULL;
 416   //these guys we don't need to copy, just reinit
 417   if (that.fQA) {fQA->Delete(); delete fQA; fQA=NULL; DefineHistograms();}
 418   fCutMC=that.fCutMC;
 419   fCutMChasTrackReferences=that.fCutMChasTrackReferences;
 420   fCutMCprocessType=that.fCutMCprocessType;
 421   fMCprocessType=that.fMCprocessType;
 422   fCutMCPID=that.fCutMCPID;
 423   fMCPID=that.fMCPID;
 424   fCutMCfirstMotherPID=that.fCutMCfirstMotherPID;
 425   fMCfirstMotherPID=that.fMCfirstMotherPID;
 426   fIgnoreSignInMCPID=that.fIgnoreSignInMCPID,
 427   fCutMCisPrimary=that.fCutMCisPrimary;
 428   fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
 429   fMCisPrimary=that.fMCisPrimary;
 430   fRequireCharge=that.fRequireCharge;
 431   fFakesAreOK=that.fFakesAreOK;
 432   fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
 433   fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
 434   fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
 435   fIgnoreTPCzRange=that.fIgnoreTPCzRange;
 436   fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
 437   fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
 438   fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
 439   fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
 440   fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
 441   fCutNClustersTPC=that.fCutNClustersTPC;
 442   fNClustersTPCMax=that.fNClustersTPCMax;
 443   fNClustersTPCMin=that.fNClustersTPCMin;
 444   fCutNClustersITS=that.fCutNClustersITS;
 445   fNClustersITSMax=that.fNClustersITSMax;
 446   fNClustersITSMin=that.fNClustersITSMin;
 447   fUseAODFilterBit=that.fUseAODFilterBit;
 448   fAODFilterBit=that.fAODFilterBit;
 449   fCutDCAToVertexXY=that.fCutDCAToVertexXY;
 450   fCutDCAToVertexZ=that.fCutDCAToVertexZ;
 451   fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
 452   fMinimalTPCdedx=that.fMinimalTPCdedx;
 453   fLinearizeVZEROresponse=that.fLinearizeVZEROresponse;
 454   fCutPmdDet=that.fCutPmdDet;
 455   fPmdDet=that.fPmdDet;
 456   fCutPmdAdc=that.fCutPmdAdc;
 457   fPmdAdc=that.fPmdAdc;
 458   fCutPmdNcell=that.fCutPmdNcell;
 459   fPmdNcell=that.fPmdNcell;
 460   fMinKinkAngle=that.fMinKinkAngle;
 461   fMinKinkRadius=that.fMinKinkRadius;
 462   fMaxKinkRadius=that.fMaxKinkRadius;
 463   fMinKinkQt=that.fMinKinkQt;
 464   fMaxKinkQt=that.fMaxKinkQt;
 465   fMinKinkInvMassKmu=that.fMinKinkInvMassKmu;
 466   fMaxKinkInvMassKmu=that.fMaxKinkInvMassKmu;
 467   fForceTPCstandalone=that.fForceTPCstandalone;
 468   fRequireKinkDaughters=that.fRequireKinkDaughters;
 469
 470   fParamType=that.fParamType;
 471   fParamMix=that.fParamMix;
 472
 473   fKink=NULL;
 474   fV0=NULL;
 475   fTrack=NULL;
 476   fTrackMass=0.;
 477   fTrackPt=0.;
 478   fTrackEta=0.;
 479   fTrackPhi=0.;
 480   fTrackWeight=1.;
 481   fTrackLabel=INT_MIN;
 482   fMCevent=NULL;
 483   fMCparticle=NULL;
 484   fEvent=NULL;
 485
 486   fESDpid = that.fESDpid;
 487   fPIDsource = that.fPIDsource;
 488
 489   delete fTPCpidCuts;
 490   delete fTOFpidCuts;
 491   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
 492   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
 493
 494   fParticleID=that.fParticleID;
 495   fParticleProbability=that.fParticleProbability;
 496   fAllowTOFmismatchFlag=that.fAllowTOFmismatchFlag;
 497   fRequireStrictTOFTPCagreement=that.fRequireStrictTOFTPCagreement;
 498   fCutRejectElectronsWithTPCpid=that.fCutRejectElectronsWithTPCpid;
 499   fProbBayes = that.fProbBayes;
 500   fCurrCentr = that.fCurrCentr;
 501
 502   fApplyRecentering = that.fApplyRecentering;
 503   fVZEROgainEqualizationPerRing = that.fVZEROgainEqualizationPerRing;
 504 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
 505   if (that.fVZEROgainEqualization) fVZEROgainEqualization = new TH1(*(that.fVZEROgainEqualization));
 506 #else
 507   //PH Lets try Clone, however the result might be wrong
 508   if (that.fVZEROgainEqualization) fVZEROgainEqualization = (TH1*)that.fVZEROgainEqualization->Clone();
 509 #endif
 510   for(Int_t i(0); i < 4; i++) { // no use to copy these guys since they're only initialized on worked node
 511       fVZEROApol[i] = that.fVZEROApol[i];
 512       fVZEROCpol[i] = that.fVZEROCpol[i];
 513   }
 514   for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = that.fUseVZERORing[i];
 515
 516   return *this;
 517 }
 518
 519 //-----------------------------------------------------------------------
 520 AliFlowTrackCuts::~AliFlowTrackCuts()
 521 {
 522   //dtor
 523   delete fAliESDtrackCuts;
 524   delete fTPCpidCuts;
 525   delete fTOFpidCuts;
 526   if (fMuonTrackCuts) delete fMuonTrackCuts;  // XZhang 20120604
 527   if (fQA) { fQA->SetOwner(); fQA->Delete(); delete fQA; }
 528   if (fVZEROgainEqualization) delete fVZEROgainEqualization;
 529 }
 530
 531 //-----------------------------------------------------------------------
 532 void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
 533 {
 534   //set the event
 535   Clear();
 536   fEvent=event;
 537   fMCevent=mcEvent;
 538
 539   //do the magic for ESD
 540   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
 541   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(event);
 542   if (fCutPID && myESD)
 543   {
 544     //TODO: maybe call it only for the TOF options?
 545     // Added by F. Noferini for TOF PID
 546     // old procedure now implemented inside fBayesianResponse
 547     //    fESDpid.MakePID(myESD,kFALSE);
 548     // new procedure
 549     fBayesianResponse->SetDetResponse(myESD, fCurrCentr,AliESDpid::kTOF_T0,kTRUE); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
 550     fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
 551     // End F. Noferini added part
 552   }
 553   if (fCutPID && myAOD){
 554     fBayesianResponse->SetDetResponse(myAOD, fCurrCentr,AliESDpid::kTOF_T0); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
 555     if(myAOD->GetTOFHeader()){
 556       fESDpid.SetTOFResponse(myAOD,AliESDpid::kTOF_T0);
 557     }
 558     else{ // corrected on the fly track by track if tof header is not present (old AODs)
 559     }
 560     // End F. Noferini added part
 561   }
 562
 563   if(fPIDsource==kTOFbayesian) fBayesianResponse->SetDetAND(1);
 564   else if(fPIDsource==kTPCbayesian) fBayesianResponse->ResetDetOR(1);
 565 }
 566
 567 //-----------------------------------------------------------------------
 568 void AliFlowTrackCuts::SetCutMC( Bool_t b )
 569 {
 570   //will we be cutting on MC information?
 571   fCutMC=b;
 572
 573   //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
 574   if (fCutMC)
 575   {
 576     fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
 577                                                        1.75295e+01,
 578                                                        3.40030e-09,
 579                                                        1.96178e+00,
 580                                                        3.91720e+00);
 581   }
 582 }
 583
 584 //-----------------------------------------------------------------------
 585 Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
 586 {
 587   //check cuts
 588   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 589 //if (vparticle) return PassesCuts(vparticle);                // XZhang 20120604
 590   if (vparticle) {                                            // XZhang 20120604
 591     if (fParamType==kMUON) return PassesMuonCuts(vparticle);  // XZhang 20120604
 592     return PassesCuts(vparticle);                             // XZhang 20120604
 593   }                                                           // XZhang 20120604
 594
 595   AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
 596   if (flowtrack) return PassesCuts(flowtrack);
 597   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 598   if (tracklets) return PassesCuts(tracklets,id);
 599   AliAODTracklets* trkletAOD = dynamic_cast<AliAODTracklets*>(obj);  // XZhang 20120615
 600   if (trkletAOD) return PassesCuts(trkletAOD,id);                    // XZhang 20120615
 601   AliESDPmdTrack* pmdtrack = dynamic_cast<AliESDPmdTrack*>(obj);
 602   if (pmdtrack) return PassesPMDcuts(pmdtrack);
 603   AliVEvent* vvzero = dynamic_cast<AliVEvent*>(obj); // should be removed; left for protection only
 604   if (vvzero) return PassesVZEROcuts(id);
 605   AliESDkink* kink = dynamic_cast<AliESDkink*>(obj);
 606   if (kink) return PassesCuts(kink);
 607   //AliESDv0* v0 = dynamic_cast<AliESDv0*>(obj);
 608   //if (v0) return PassesCuts(v0);
 609
 610   return kFALSE;  //default when passed wrong type of object
 611 }
 612
 613 //-----------------------------------------------------------------------
 614 Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
 615 {
 616   //check cuts
 617   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
 618   if (vparticle)
 619   {
 620     return PassesMCcuts(fMCevent,(fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel());
 621   }
 622   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
 623   if (tracklets)
 624   {
 625     Int_t label0 = tracklets->GetLabel(id,0);
 626     Int_t label1 = tracklets->GetLabel(id,1);
 627     Int_t label = (fFakesAreOK)?TMath::Abs(label0):label0;
 628     if (label0!=label1) label = -666;
 629     return PassesMCcuts(fMCevent,label);
 630   }
 631   return kFALSE;  //default when passed wrong type of object
 632 }
 633
 634 //-----------------------------------------------------------------------
 635 Bool_t AliFlowTrackCuts::PassesCuts(const AliFlowTrackSimple* track)
 636 {
 637   //check cuts on a flowtracksimple
 638
 639   //clean up from last iteration
 640   ClearTrack();
 641   return AliFlowTrackSimpleCuts::PassesCuts(track);
 642 }
 643
 644 //-----------------------------------------------------------------------
 645 Bool_t AliFlowTrackCuts::PassesCuts(const AliMultiplicity* tracklet, Int_t id)
 646 {
 647   //check cuts on a tracklets
 648
 649   if (id<0) return kFALSE;
 650
 651   //clean up from last iteration, and init label
 652   ClearTrack();
 653
 654   fTrackPhi = tracklet->GetPhi(id);
 655   fTrackEta = tracklet->GetEta(id);
 656   fTrackWeight = 1.0;
 657   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
 658   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
 659
 660   //check MC info if available
 661   //if the 2 clusters have different label track cannot be good
 662   //and should therefore not pass the mc cuts
 663   Int_t label0 = tracklet->GetLabel(id,0);
 664   Int_t label1 = tracklet->GetLabel(id,1);
 665   //if possible get label and mcparticle
 666   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
 667   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
 668   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 669   //check MC cuts
 670   if (fCutMC && !PassesMCcuts()) return kFALSE;
 671   return kTRUE;
 672 }
 673
 674 //-----------------------------------------------------------------------
 675 Bool_t AliFlowTrackCuts::PassesCuts(const AliAODTracklets* tracklet, Int_t id)
 676 {
 677   // XZhang 20120615
 678   //check cuts on a tracklets
 679
 680   if (id<0) return kFALSE;
 681
 682   //clean up from last iteration, and init label
 683   ClearTrack();
 684
 685   fTrackPhi = tracklet->GetPhi(id);
 686 //fTrackEta = tracklet->GetEta(id);
 687   fTrackEta = -1.*TMath::Log(TMath::Tan(tracklet->GetTheta(id)/2.));
 688   fTrackWeight = 1.0;
 689   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
 690   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
 691
 692   //check MC info if available
 693   //if the 2 clusters have different label track cannot be good
 694   //and should therefore not pass the mc cuts
 695   Int_t label0 = tracklet->GetLabel(id,0);
 696   Int_t label1 = tracklet->GetLabel(id,1);
 697   //if possible get label and mcparticle
 698   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
 699   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
 700   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 701   //check MC cuts
 702   if (fCutMC && !PassesMCcuts()) return kFALSE;
 703   return kTRUE;
 704 }
 705
 706 //-----------------------------------------------------------------------
 707 Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
 708 {
 709   //check the MC info
 710   if (!mcEvent) return kFALSE;
 711   if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 712   AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
 713   if (!mcparticle) {AliError("no MC track"); return kFALSE;}
 714
 715   if (fCutMCisPrimary)
 716   {
 717     if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
 718   }
 719   if (fCutMCPID)
 720   {
 721     Int_t pdgCode = mcparticle->PdgCode();
 722     if (fIgnoreSignInMCPID)
 723     {
 724       if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
 725     }
 726     else
 727     {
 728       if (fMCPID != pdgCode) return kFALSE;
 729     }
 730   }
 731   if (fCutMCfirstMotherPID)
 732   {
 733
 734     TParticle* tparticle=mcparticle->Particle();
 735     Int_t firstMotherLabel = 0;
 736     if (tparticle) { firstMotherLabel = tparticle->GetFirstMother(); }
 737     AliVParticle* firstMotherParticle = mcEvent->GetTrack(firstMotherLabel);
 738     Int_t pdgcodeFirstMother = 0;
 739     if (firstMotherParticle) { pdgcodeFirstMother = firstMotherParticle->PdgCode(); }
 740     if (pdgcodeFirstMother != fMCfirstMotherPID) return kFALSE;
 741   }
 742   if ( fCutMCprocessType )
 743   {
 744     TParticle* particle = mcparticle->Particle();
 745     Int_t processID = particle->GetUniqueID();
 746     if (processID != fMCprocessType ) return kFALSE;
 747   }
 748   if (fCutMChasTrackReferences)
 749   {
 750     if (mcparticle->GetNumberOfTrackReferences()<1) return kFALSE;
 751   }
 752   return kTRUE;
 753 }
 754
 755 //-----------------------------------------------------------------------
 756 Bool_t AliFlowTrackCuts::PassesMCcuts()
 757 {
 758   //check MC info
 759   if (!fMCevent) return kFALSE;
 760   if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
 761   fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 762   return PassesMCcuts(fMCevent,fTrackLabel);
 763 }
 764
 765 //-----------------------------------------------------------------------
 766 Bool_t AliFlowTrackCuts::PassesCuts(const AliESDv0* v0)
 767 {
 768   //check if the v0 passes cuts
 769
 770   //clean up from last iteration
 771   ClearTrack();
 772
 773   fV0 = const_cast<AliESDv0*>(v0);
 774
 775   Bool_t pass = kTRUE;
 776
 777   //first, extract/prepare the v0
 778   if (!v0->GetOnFlyStatus()) return kFALSE; //skip offline V0 finder
 779   const AliExternalTrackParam *negHelix=v0->GetParamN();
 780   const AliExternalTrackParam *posHelix=v0->GetParamP();
 781   AliVParticle *v0tracks[2];
 782   v0tracks[0] = fEvent->GetTrack(v0->GetNindex());
 783   v0tracks[1] = fEvent->GetTrack(v0->GetPindex());
 784   if( v0tracks[1]->Charge() < 0)
 785   {
 786     v0tracks[1] = fEvent->GetTrack(v0->GetNindex());
 787     v0tracks[0] = fEvent->GetTrack(v0->GetPindex());
 788     negHelix=v0->GetParamP();
 789     posHelix=v0->GetParamN();
 790   }
 791
 792   int KalmanPidPairs[4][2] =
 793   {
 794     {-11,11}, // e-,e+ (gamma)
 795     {-211,2212}, // pi- p (lambda)
 796     {-2212,211}, // anti-p pi+ (anti-lambda)
 797     {-211,211} // pi-  pi+ (Kshort)
 798     //   {-321,321} // K- K+ (phi)
 799   };
 800
 801   Int_t id = 3;
 802   //refit using a mass hypothesis
 803   AliKFParticle v0trackKFneg(*(negHelix),KalmanPidPairs[id][0]);
 804   AliKFParticle v0trackKFpos(*(posHelix),KalmanPidPairs[id][1]);
 805   AliKFParticle v0particleRefit;
 806   v0particleRefit += v0trackKFneg;
 807   v0particleRefit += v0trackKFpos;
 808   Double_t invMassErr= -999;
 809   v0particleRefit.GetMass(fTrackMass,invMassErr);
 810   //Double_t openAngle = v0trackKFneg.GetAngle(v0trackKFpos);
 811   fTrackEta = v0particleRefit.GetEta();
 812   fTrackPt = v0particleRefit.GetPt();
 813   fTrackPhi = TMath::Pi()+v0particleRefit.GetPhi();
 814   ////find out which mass bin and put the number in fPOItype
 815   //Int_t massBins = AliFlowCommonConstants::GetMaster()->GetNbinsMass();
 816   //Double_t massMin = AliFlowCommonConstants::GetMaster()->GetMassMin();
 817   //Double_t massMax = AliFlowCommonConstants::GetMaster()->GetMassMax();
 818   //fPOItype = 1 + int(massBins*(fTrackMass-massMin)/(massMax-massMin));
 819
 820   /////////////////////////////////////////////////////////////////////////////
 821   //apply cuts
 822   if ( v0tracks[0]->Charge() == v0tracks[1]->Charge() ) pass=kFALSE;
 823   if ( v0tracks[0]->Pt()<0.15 || v0tracks[1]->Pt()<0.15 ) pass=kFALSE;
 824
 825   return pass;
 826 }
 827
 828 //-----------------------------------------------------------------------
 829 Bool_t AliFlowTrackCuts::PassesCuts(const AliESDkink* kink)
 830 {
 831   //check if the kink passes cuts
 832
 833   //clean up from last iteration
 834   ClearTrack();
 835   fKink=const_cast<AliESDkink*>(kink);
 836
 837   Bool_t pass = kTRUE;
 838
 839   Float_t kinkAngle = kink->GetAngle(2);
 840   if (kinkAngle<fMinKinkAngle) pass = kFALSE;
 841   Double_t kinkRadius = kink->GetR();
 842   if (kinkRadius<fMinKinkRadius || kinkRadius>fMaxKinkRadius) pass = kFALSE;
 843
 844   //Qt
 845   const TVector3 motherMfromKink(kink->GetMotherP());
 846   const TVector3 daughterMfromKink(kink->GetDaughterP());
 847   Float_t qt=kink->GetQt();
 848   if ( qt < fMinKinkQt || qt > fMaxKinkQt) pass = kFALSE;
 849
 850   //invariant mass
 851   Float_t energyDaughterMu = TMath::Sqrt( daughterMfromKink.Mag()*daughterMfromKink.Mag()+
 852                                           0.105658*0.105658 );
 853   Float_t p1XM = motherMfromKink.Px();
 854   Float_t p1YM = motherMfromKink.Py();
 855   Float_t p1ZM = motherMfromKink.Pz();
 856   Float_t p2XM = daughterMfromKink.Px();
 857   Float_t p2YM = daughterMfromKink.Py();
 858   Float_t p2ZM = daughterMfromKink.Pz();
 859   Float_t p3Daughter = TMath::Sqrt( ((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+
 860                                     ((p1ZM-p2ZM)*(p1ZM-p2ZM)) );
 861   Double_t invariantMassKmu = TMath::Sqrt( (energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-
 862                                            motherMfromKink.Mag()*motherMfromKink.Mag() );
 863   if (invariantMassKmu > fMaxKinkInvMassKmu) pass=kFALSE;
 864   if (invariantMassKmu < fMinKinkInvMassKmu) pass=kFALSE;
 865   fTrackMass=invariantMassKmu;
 866
 867   if (fQA)
 868   {
 869     QAbefore(13)->Fill(qt);
 870     if (pass) QAafter(13)->Fill(qt);
 871     QAbefore(14)->Fill(invariantMassKmu);
 872     if (pass) QAafter(14)->Fill(invariantMassKmu);
 873     const Double_t* kinkPosition = kink->GetPosition();
 874     QAbefore(15)->Fill(kinkPosition[0],kinkPosition[1]);
 875     if (pass) QAafter(15)->Fill(kinkPosition[0],kinkPosition[1]);
 876     QAbefore(16)->Fill(motherMfromKink.Mag(),kinkAngle*TMath::RadToDeg());
 877     if (pass) QAafter(16)->Fill(motherMfromKink.Mag(),kinkAngle*TMath::RadToDeg());
 878   }
 879
 880   //mother track cuts
 881   Int_t indexKinkMother = kink->GetIndex(0);
 882   AliESDtrack* motherTrack = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(indexKinkMother));
 883   if (!motherTrack) return kFALSE;
 884   if (!PassesCuts(motherTrack)) pass = kFALSE;
 885
 886   return pass;
 887 }
 888
 889 //-----------------------------------------------------------------------
 890 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
 891 {
 892   //check cuts for an ESD vparticle
 893
 894   ////////////////////////////////////////////////////////////////
 895   //  start by preparing the track parameters to cut on //////////
 896   ////////////////////////////////////////////////////////////////
 897   //clean up from last iteration
 898   ClearTrack();
 899   Bool_t pass=kTRUE;
 900
 901   //get the label and the mc particle
 902   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
 903   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
 904   else fMCparticle=NULL;
 905
 906   Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
 907   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
 908   AliAODTrack* aodTrack = NULL;
 909   if (esdTrack)
 910   {
 911     //for an ESD track we do some magic sometimes like constructing TPC only parameters
 912     //or doing some hybrid, handle that here
 913     HandleESDtrack(esdTrack);
 914   }
 915   else
 916   {
 917     HandleVParticle(vparticle);
 918     //now check if produced particle is MC
 919     isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
 920     aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
 921   }
 922   ////////////////////////////////////////////////////////////////
 923   ////////////////////////////////////////////////////////////////
 924
 925   if (!fTrack) return kFALSE;
 926   //because it may be different from global, not needed for aodTrack because we dont do anything funky there
 927   if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
 928
 929   //check the common cuts for the current particle fTrack (MC,AOD,ESD)
 930   Double_t pt = fTrack->Pt();
 931   Double_t p = fTrack->P();
 932   if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
 933   if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
 934   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
 935   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
 936   if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
 937   if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
 938   if (fCutCharge && isMCparticle)
 939   {
 940     //in case of an MC particle the charge is stored in units of 1/3|e|
 941     Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
 942     if (charge!=fCharge) pass=kFALSE;
 943   }
 944   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
 945
 946   //when additionally MC info is required
 947   if (fCutMC && !PassesMCcuts()) pass=kFALSE;
 948
 949   //the case of ESD or AOD
 950   if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
 951   if (aodTrack) { if (!PassesAODcuts(aodTrack,pass)) { pass=kFALSE; } }
 952
 953   if (fQA)
 954   {
 955     if (fMCparticle)
 956     {
 957       TParticle* tparticle=fMCparticle->Particle();
 958       Int_t processID = tparticle->GetUniqueID();
 959       Int_t firstMotherLabel = tparticle->GetFirstMother();
 960       Bool_t primary = IsPhysicalPrimary();
 961       //TLorentzVector v;
 962       //mcparticle->Particle()->ProductionVertex(v);
 963       //Double_t prodvtxX = v.X();
 964       //Double_t prodvtxY = v.Y();
 965
 966       Int_t pdgcode = fMCparticle->PdgCode();
 967       AliVParticle* firstMotherParticle = fMCevent->GetTrack(firstMotherLabel);
 968       Int_t pdgcodeFirstMother = 0;
 969       if (firstMotherParticle) {pdgcodeFirstMother = firstMotherParticle->PdgCode();}
 970
 971       Float_t pdg = 0;
 972       switch (TMath::Abs(pdgcode))
 973       {
 974         case 11:
 975           pdg = AliPID::kElectron + 0.5; break;
 976         case 13:
 977           pdg = AliPID::kMuon + 0.5; break;
 978         case 211:
 979           pdg = AliPID::kPion + 0.5; break;
 980         case 321:
 981           pdg = AliPID::kKaon + 0.5; break;
 982         case 2212:
 983           pdg = AliPID::kProton + 0.5; break;
 984         default:
 985           pdg = AliPID::kUnknown + 0.5; break;
 986       }
 987       pdg = TMath::Sign(pdg,static_cast<Float_t>(pdgcode));
 988
 989       Float_t geantCode = 0;
 990       switch (pdgcodeFirstMother)
 991       {
 992         case 22: //#gamma
 993           geantCode=1;
 994           break;
 995         case -11: //e^{+}
 996           geantCode=2;
 997           break;
 998         case 11: //e^{-}
 999           geantCode=3;
1000           break;
1001         case 12: case 14: case 16: //#nu
1002           geantCode=4;
1003           break;
1004         case -13:
1005           geantCode=5; //#mu^{+}
1006           break;
1007         case 13:
1008           geantCode=6; //#mu^{-}
1009           break;
1010         case 111: //#pi^{0}
1011           geantCode=7;
1012           break;
1013         case 211: //#pi^{+}
1014           geantCode=8;
1015           break;
1016         case -211: //#pi^{-}
1017           geantCode=9;
1018           break;
1019         case 130: //K^{0}_{L}
1020           geantCode=10;
1021           break;
1022         case 321: //K^{+}
1023           geantCode=11;
1024           break;
1025         case -321: //K^{-}
1026           geantCode=12;
1027           break;
1028         case 2112:
1029           geantCode = 13; //n
1030           break;
1031         case 2212:
1032           geantCode = 14; //p
1033           break;
1034         case -2212:
1035           geantCode=15; //#bar{p}
1036           break;
1037         case 310:
1038           geantCode=16; //K^{0}_{S}
1039           break;
1040         case 221:
1041           geantCode=17; //#eta
1042           break;
1043         case 3122:
1044           geantCode=18; //#Lambda
1045           break;
1046         case 3222:
1047           geantCode=19; //#Sigma^{+}
1048           break;
1049         case 3212:
1050           geantCode=20; //#Sigma^{0}
1051           break;
1052         case 3112:
1053           geantCode=21; //#Sigma^{-}
1054           break;
1055         case 3322:
1056           geantCode=22; //#Theta^{0}
1057           break;
1058         case 3312:
1059           geantCode=23; //#Theta^{-}
1060           break;
1061         case 3332:
1062           geantCode=24; //#Omega^{-}
1063           break;
1064         case -2112:
1065           geantCode=25; //#bar{n}
1066           break;
1067         case -3122:
1068           geantCode=26; //#bar{#Lambda}
1069           break;
1070         case -3112:
1071           geantCode=27; //#bar{#Sigma}^{-}
1072           break;
1073         case -3212:
1074           geantCode=28; //#bar{#Sigma}^{0}
1075           break;
1076         case -3222:
1077           geantCode=29; //#bar{#Sigma}^{+}
1078           break;
1079         case -3322:
1080           geantCode=30; //#bar{#Theta}^{0}
1081           break;
1082         case -3312:
1083           geantCode=31; //#bar{#Theta}^{+}
1084           break;
1085         case -3332:
1086           geantCode=32; //#bar{#Omega}^{+}
1087           break;
1088         case -15:
1089           geantCode=33; //#tau^{+}
1090           break;
1091         case 15:
1092           geantCode=34; //#tau^{-}
1093           break;
1094         case 411:
1095           geantCode=35; //D^{+}
1096           break;
1097         case -411:
1098           geantCode=36; //D^{-}
1099           break;
1100         case 421:
1101           geantCode=37; //D^{0}
1102           break;
1103         case -421:
1104           geantCode=38; //#bar{D}^{0}
1105           break;
1106         case 431:
1107           geantCode=39; //D_{s}^{+}
1108           break;
1109         case -431:
1110           geantCode=40; //#bar{D_{s}}^{-}
1111           break;
1112         case 4122:
1113           geantCode=41; //#Lambda_{c}^{+}
1114           break;
1115         case 24:
1116           geantCode=42; //W^{+}
1117           break;
1118         case -24:
1119           geantCode=43; //W^{-}
1120           break;
1121         case 23:
1122           geantCode=44; //Z^{0}
1123           break;
1124         default:
1125           geantCode = 100;
1126           break;
1127       }
1128
1129       QAbefore(2)->Fill(p,pdg);
1130       QAbefore(3)->Fill(p,primary?0.5:-0.5);
1131       QAbefore(4)->Fill(p,static_cast<Float_t>(processID));
1132       QAbefore(7)->Fill(p,geantCode+0.5);
1133       if (pass) QAafter(2)->Fill(p,pdg);
1134       if (pass) QAafter(3)->Fill(p,primary?0.5:-0.5);
1135       if (pass) QAafter(4)->Fill(p,static_cast<Float_t>(processID));
1136       if (pass) QAafter(7)->Fill(p,geantCode);
1137     }
1138   }
1139
1140   //true by default, if we didn't set any cuts
1141   return pass;
1142 }
1143
1144 //_______________________________________________________________________
1145 Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track, Bool_t passedFid)
1146 {
1147   //check cuts for AOD
1148   Bool_t pass = passedFid;
1149
1150   if (fCutNClustersTPC)
1151   {
1152     Int_t ntpccls = track->GetTPCNcls();
1153     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1154   }
1155
1156   if (fCutNClustersITS)
1157   {
1158     Int_t nitscls = track->GetITSNcls();
1159     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1160   }
1161
1162    if (fCutChi2PerClusterTPC)
1163   {
1164     Double_t chi2tpc = track->Chi2perNDF();
1165     if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
1166   }
1167
1168   if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
1169   if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
1170
1171   if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
1172
1173   if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
1174
1175   if (fCutDCAToVertexZ && track->ZAtDCA()>GetMaxDCAToVertexZ()) pass=kFALSE;
1176
1177   Double_t dedx = track->GetTPCsignal();
1178   if (dedx < fMinimalTPCdedx) pass=kFALSE;
1179   Double_t time[9];
1180   track->GetIntegratedTimes(time);
1181   if (fQA) {
1182     Double_t momTPC = track->GetTPCmomentum();
1183     QAbefore( 1)->Fill(momTPC,dedx);
1184     QAbefore( 5)->Fill(track->Pt(),track->DCA());
1185     QAbefore( 6)->Fill(track->Pt(),track->ZAtDCA());
1186     if (pass) QAafter( 1)->Fill(momTPC,dedx);
1187     if (pass) QAafter( 5)->Fill(track->Pt(),track->DCA());
1188     if (pass) QAafter( 6)->Fill(track->Pt(),track->ZAtDCA());
1189     QAbefore( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1190     if (pass) QAafter(  8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1191     QAbefore( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1192     if (pass) QAafter(  9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1193     QAbefore(10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1194     if (pass) QAafter( 10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1195     QAbefore(11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1196     if (pass) QAafter( 11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1197     QAbefore(12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1198     if (pass) QAafter( 12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1199   }
1200
1201   if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1202     {
1203       if (!PassesAODpidCut(track)) pass=kFALSE;
1204     }
1205
1206   return pass;
1207 }
1208
1209 //_______________________________________________________________________
1210 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
1211 {
1212   //check cuts on ESD tracks
1213   Bool_t pass=kTRUE;
1214   Float_t dcaxy=0.0;
1215   Float_t dcaz=0.0;
1216   track->GetImpactParameters(dcaxy,dcaz);
1217   const AliExternalTrackParam* pout = track->GetOuterParam();
1218   const AliExternalTrackParam* pin = track->GetInnerParam();
1219   if (fIgnoreTPCzRange)
1220   {
1221     if (pin&&pout)
1222     {
1223       Double_t zin = pin->GetZ();
1224       Double_t zout = pout->GetZ();
1225       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
1226       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
1227       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
1228     }
1229   }
1230
1231   Int_t ntpccls = ( fParamType==kTPCstandalone )?
1232                     track->GetTPCNclsIter1():track->GetTPCNcls();
1233   if (fCutChi2PerClusterTPC)
1234   {
1235     Float_t tpcchi2 = (fParamType==kTPCstandalone)?
1236                        track->GetTPCchi2Iter1():track->GetTPCchi2();
1237     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
1238     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
1239       pass=kFALSE;
1240   }
1241
1242   if (fCutMinimalTPCdedx)
1243   {
1244     if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
1245   }
1246
1247   if (fCutNClustersTPC)
1248   {
1249     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1250   }
1251
1252   Int_t nitscls = track->GetNcls(0);
1253   if (fCutNClustersITS)
1254   {
1255     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1256   }
1257
1258   //some stuff is still handled by AliESDtrackCuts class - delegate
1259   if (fAliESDtrackCuts)
1260   {
1261     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
1262   }
1263
1264   //PID part with pid QA
1265   Double_t beta = GetBeta(track);
1266   Double_t dedx = Getdedx(track);
1267   if (fQA)
1268   {
1269     if (pass) QAbefore(0)->Fill(track->GetP(),beta);
1270     if (pass) QAbefore(1)->Fill(pin->GetP(),dedx);
1271   }
1272   if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1273   {
1274     if (!PassesESDpidCut(track)) pass=kFALSE;
1275   }
1276   if (fCutRejectElectronsWithTPCpid)
1277   {
1278     //reject electrons using standard TPC pid
1279     //TODO this should be rethought....
1280     Double_t pidTPC[AliPID::kSPECIES];
1281     track->GetTPCpid(pidTPC);
1282     if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
1283   }
1284   if (fQA)
1285   {
1286     if (pass) QAafter(0)->Fill(track->GetP(),beta);
1287     if (pass) QAafter(1)->Fill(pin->GetP(),dedx);
1288   }
1289   //end pid part with qa
1290
1291   if (fQA)
1292   {
1293     Double_t pt = track->Pt();
1294     QAbefore(5)->Fill(pt,dcaxy);
1295     QAbefore(6)->Fill(pt,dcaz);
1296     if (pass) QAafter(5)->Fill(pt,dcaxy);
1297     if (pass) QAafter(6)->Fill(pt,dcaz);
1298     QAbefore(17)->Fill(Float_t(track->GetKinkIndex(0)));
1299     if (pass) QAafter(17)->Fill(Float_t(track->GetKinkIndex(0)));
1300   }
1301
1302   return pass;
1303 }
1304
1305 //-----------------------------------------------------------------------
1306 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
1307 {
1308   //handle the general case
1309   switch (fParamType)
1310   {
1311     default:
1312       fTrack = track;
1313       break;
1314   }
1315 }
1316
1317 //-----------------------------------------------------------------------
1318 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
1319 {
1320   //handle esd track
1321   switch (fParamType)
1322   {
1323     case kGlobal:
1324       fTrack = track;
1325       break;
1326     case kTPCstandalone:
1327       if (!track->FillTPCOnlyTrack(fTPCtrack))
1328       {
1329         fTrack=NULL;
1330         fMCparticle=NULL;
1331         fTrackLabel=-998;
1332         return;
1333       }
1334       fTrack = &fTPCtrack;
1335       //recalculate the label and mc particle, they may differ as TPClabel != global label
1336       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1337       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1338       else fMCparticle=NULL;
1339       break;
1340     default:
1341       if (fForceTPCstandalone)
1342       {
1343         if (!track->FillTPCOnlyTrack(fTPCtrack))
1344         {
1345           fTrack=NULL;
1346           fMCparticle=NULL;
1347           fTrackLabel=-998;
1348           return;
1349         }
1350         fTrack = &fTPCtrack;
1351         //recalculate the label and mc particle, they may differ as TPClabel != global label
1352         fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1353         if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1354         else fMCparticle=NULL;
1355       }
1356       else
1357         fTrack = track;
1358       break;
1359   }
1360 }
1361
1362 //-----------------------------------------------------------------------
1363 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
1364 {
1365   //calculate the number of track in given event.
1366   //if argument is NULL(default) take the event attached
1367   //by SetEvent()
1368   Int_t multiplicity = 0;
1369   if (!event)
1370   {
1371     for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
1372     {
1373       if (IsSelected(GetInputObject(i))) multiplicity++;
1374     }
1375   }
1376   else
1377   {
1378     for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
1379     {
1380       if (IsSelected(event->GetTrack(i))) multiplicity++;
1381     }
1382   }
1383   return multiplicity;
1384 }
1385
1386 //-----------------------------------------------------------------------
1387 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts()
1388 {
1389   //returns the lhc10h vzero track cuts, this function
1390   //is left here for backward compatibility
1391   //if a run is recognized as 11h, the calibration method will
1392   //switch to 11h calbiration, which means that the cut
1393   //object is updated but not replaced.
1394   //calibratin is only available for PbPb runs
1395   return GetStandardVZEROOnlyTrackCuts2010();
1396 }
1397 //-----------------------------------------------------------------------
1398 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010()
1399 {
1400   //get standard VZERO cuts
1401   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2010");
1402   cuts->SetParamType(kVZERO);
1403   cuts->SetEtaRange( -10, +10 );
1404   cuts->SetPhiMin( 0 );
1405   cuts->SetPhiMax( TMath::TwoPi() );
1406   // options for the reweighting
1407   cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1408   cuts->SetApplyRecentering(kTRUE);
1409   // to exclude a ring , do e.g.
1410   // cuts->SetUseVZERORing(7, kFALSE);
1411   return cuts;
1412 }
1413 //-----------------------------------------------------------------------
1414 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011()
1415 {
1416   //get standard VZERO cuts for 2011 data
1417   //in this case, the vzero segments will be weighted by
1418   //VZEROEqMultiplicity,
1419   //if recentering is enableded, the sub-q vectors
1420   //will be taken from the event header, so make sure to run
1421   //the VZERO event plane selection task before this task !
1422   //recentering replaces the already evaluated q-vectors, so
1423   //when chosen, additional settings (e.g. excluding rings)
1424   //have no effect. recentering is true by default
1425   //
1426   //NOTE user is responsible for running the vzero event plane
1427   //selection task in advance, e.g. add to your launcher macro
1428   //
1429   //  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C");
1430   //  AddTaskVZEROEPSelection();
1431   //
1432   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2011");
1433   cuts->SetParamType(kVZERO);
1434   cuts->SetEtaRange( -10, +10 );
1435   cuts->SetPhiMin( 0 );
1436   cuts->SetPhiMax( TMath::TwoPi() );
1437   cuts->SetApplyRecentering(kTRUE);
1438   return cuts;
1439 }
1440 //-----------------------------------------------------------------------
1441 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
1442 {
1443   //get standard cuts
1444   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global tracks");
1445   cuts->SetParamType(kGlobal);
1446   cuts->SetPtRange(0.2,5.);
1447   cuts->SetEtaRange(-0.8,0.8);
1448   cuts->SetMinNClustersTPC(70);
1449   cuts->SetMinChi2PerClusterTPC(0.1);
1450   cuts->SetMaxChi2PerClusterTPC(4.0);
1451   cuts->SetMinNClustersITS(2);
1452   cuts->SetRequireITSRefit(kTRUE);
1453   cuts->SetRequireTPCRefit(kTRUE);
1454   cuts->SetMaxDCAToVertexXY(0.3);
1455   cuts->SetMaxDCAToVertexZ(0.3);
1456   cuts->SetAcceptKinkDaughters(kFALSE);
1457   cuts->SetMinimalTPCdedx(10.);
1458
1459   return cuts;
1460 }
1461
1462 //-----------------------------------------------------------------------
1463 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()
1464 {
1465   //get standard cuts
1466   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone 2010");
1467   cuts->SetParamType(kTPCstandalone);
1468   cuts->SetPtRange(0.2,5.);
1469   cuts->SetEtaRange(-0.8,0.8);
1470   cuts->SetMinNClustersTPC(70);
1471   cuts->SetMinChi2PerClusterTPC(0.2);
1472   cuts->SetMaxChi2PerClusterTPC(4.0);
1473   cuts->SetMaxDCAToVertexXY(3.0);
1474   cuts->SetMaxDCAToVertexZ(3.0);
1475   cuts->SetDCAToVertex2D(kTRUE);
1476   cuts->SetAcceptKinkDaughters(kFALSE);
1477   cuts->SetMinimalTPCdedx(10.);
1478
1479   return cuts;
1480 }
1481
1482 //-----------------------------------------------------------------------
1483 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts()
1484 {
1485   //get standard cuts
1486   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone");
1487   cuts->SetParamType(kTPCstandalone);
1488   cuts->SetPtRange(0.2,5.);
1489   cuts->SetEtaRange(-0.8,0.8);
1490   cuts->SetMinNClustersTPC(70);
1491   cuts->SetMinChi2PerClusterTPC(0.2);
1492   cuts->SetMaxChi2PerClusterTPC(4.0);
1493   cuts->SetMaxDCAToVertexXY(3.0);
1494   cuts->SetMaxDCAToVertexZ(3.0);
1495   cuts->SetDCAToVertex2D(kTRUE);
1496   cuts->SetAcceptKinkDaughters(kFALSE);
1497   cuts->SetMinimalTPCdedx(10.);
1498
1499   return cuts;
1500 }
1501
1502 //-----------------------------------------------------------------------
1503 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
1504 {
1505   //get standard cuts
1506   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
1507   delete cuts->fAliESDtrackCuts;
1508   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
1509   cuts->SetParamType(kGlobal);
1510   return cuts;
1511 }
1512
1513 //-----------------------------------------------------------------------------
1514 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardMuonTrackCuts(Bool_t isMC, Int_t passN)
1515 {
1516 // XZhang 20120604
1517   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard muon track cuts");
1518   cuts->SetParamType(kMUON);
1519   cuts->SetStandardMuonTrackCuts();
1520   cuts->SetIsMuonMC(isMC);
1521   cuts->SetMuonPassNumber(passN);
1522   return cuts;
1523 }
1524
1525 //-----------------------------------------------------------------------
1526 //AliFlowTrack* AliFlowTrackCuts::FillFlowTrackV0(TObjArray* trackCollection, Int_t trackIndex) const
1527 //{
1528 //  //fill flow track from a reconstructed V0 (topological)
1529 //  AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1530 //  if (flowtrack)
1531 //  {
1532 //    flowtrack->Clear();
1533 //  }
1534 //  else
1535 //  {
1536 //    flowtrack = new AliFlowCandidateTrack();
1537 //    trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1538 //  }
1539 //
1540 //  TParticle *tmpTParticle=NULL;
1541 //  AliMCParticle* tmpAliMCParticle=NULL;
1542 //  AliExternalTrackParam* externalParams=NULL;
1543 //  AliESDtrack* esdtrack=NULL;
1544 //  switch(fParamMix)
1545 //  {
1546 //    case kPure:
1547 //      flowtrack->Set(fTrack);
1548 //      break;
1549 //    case kTrackWithMCkine:
1550 //      flowtrack->Set(fMCparticle);
1551 //      break;
1552 //    case kTrackWithMCPID:
1553 //      flowtrack->Set(fTrack);
1554 //      //flowtrack->setPID(...) from mc, when implemented
1555 //      break;
1556 //    case kTrackWithMCpt:
1557 //      if (!fMCparticle) return NULL;
1558 //      flowtrack->Set(fTrack);
1559 //      flowtrack->SetPt(fMCparticle->Pt());
1560 //      break;
1561 //    case kTrackWithPtFromFirstMother:
1562 //      if (!fMCparticle) return NULL;
1563 //      flowtrack->Set(fTrack);
1564 //      tmpTParticle = fMCparticle->Particle();
1565 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1566 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
1567 //      break;
1568 //    case kTrackWithTPCInnerParams:
1569 //      esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1570 //      if (!esdtrack) return NULL;
1571 //      externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1572 //      if (!externalParams) return NULL;
1573 //      flowtrack->Set(externalParams);
1574 //      break;
1575 //    default:
1576 //      flowtrack->Set(fTrack);
1577 //      break;
1578 //  }
1579 //  if (fParamType==kMC)
1580 //  {
1581 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
1582 //    flowtrack->SetID(fTrackLabel);
1583 //  }
1584 //  else if (dynamic_cast<AliAODTrack*>(fTrack))
1585 //  {
1586 //    if (fParamType==kMUON)                            // XZhang 20120604
1587 //      flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1588 //    else                                              // XZhang 20120604
1589 //      flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1590 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1591 //  }
1592 //  else if (dynamic_cast<AliMCParticle*>(fTrack))
1593 //  {
1594 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
1595 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1596 //  }
1597 //
1598 //  if (fV0)
1599 //  {
1600 //    //add daughter indices
1601 //  }
1602 //
1603 //  return flowtrack;
1604 //}
1605
1606 //-----------------------------------------------------------------------
1607 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackKink(TObjArray* trackCollection, Int_t trackIndex) const
1608 {
1609   //fill flow track from AliVParticle (ESD,AOD,MC)
1610   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1611   if (flowtrack)
1612   {
1613     flowtrack->Clear();
1614   }
1615   else
1616   {
1617     flowtrack = new AliFlowCandidateTrack();
1618     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1619   }
1620
1621   TParticle *tmpTParticle=NULL;
1622   AliMCParticle* tmpAliMCParticle=NULL;
1623   AliExternalTrackParam* externalParams=NULL;
1624   AliESDtrack* esdtrack=NULL;
1625   switch(fParamMix)
1626   {
1627     case kPure:
1628       flowtrack->Set(fTrack);
1629       break;
1630     case kTrackWithMCkine:
1631       flowtrack->Set(fMCparticle);
1632       break;
1633     case kTrackWithMCPID:
1634       flowtrack->Set(fTrack);
1635       //flowtrack->setPID(...) from mc, when implemented
1636       break;
1637     case kTrackWithMCpt:
1638       if (!fMCparticle) return NULL;
1639       flowtrack->Set(fTrack);
1640       flowtrack->SetPt(fMCparticle->Pt());
1641       break;
1642     case kTrackWithPtFromFirstMother:
1643       if (!fMCparticle) return NULL;
1644       flowtrack->Set(fTrack);
1645       tmpTParticle = fMCparticle->Particle();
1646       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1647       flowtrack->SetPt(tmpAliMCParticle->Pt());
1648       break;
1649     case kTrackWithTPCInnerParams:
1650       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1651       if (!esdtrack) return NULL;
1652       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1653       if (!externalParams) return NULL;
1654       flowtrack->Set(externalParams);
1655       break;
1656     default:
1657       flowtrack->Set(fTrack);
1658       break;
1659   }
1660   if (fParamType==kMC)
1661   {
1662     flowtrack->SetSource(AliFlowTrack::kFromMC);
1663     flowtrack->SetID(fTrackLabel);
1664   }
1665   else if (dynamic_cast<AliESDtrack*>(fTrack))
1666   {
1667     flowtrack->SetSource(AliFlowTrack::kFromESD);
1668     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1669   }
1670   else if (dynamic_cast<AliESDMuonTrack*>(fTrack))                                  // XZhang 20120604
1671   {                                                                                 // XZhang 20120604
1672     flowtrack->SetSource(AliFlowTrack::kFromMUON);                                  // XZhang 20120604
1673     flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID());  // XZhang 20120604
1674   }                                                                                 // XZhang 20120604
1675   else if (dynamic_cast<AliAODTrack*>(fTrack))
1676   {
1677     if (fParamType==kMUON)                            // XZhang 20120604
1678       flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1679     else                                              // XZhang 20120604
1680       flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1681     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1682   }
1683   else if (dynamic_cast<AliMCParticle*>(fTrack))
1684   {
1685     flowtrack->SetSource(AliFlowTrack::kFromMC);
1686     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1687   }
1688
1689   if (fKink)
1690   {
1691     Int_t indexMother = fKink->GetIndex(0);
1692     Int_t indexDaughter = fKink->GetIndex(1);
1693     flowtrack->SetID(indexMother);
1694     flowtrack->AddDaughter(indexDaughter);
1695     flowtrack->SetMass(1.);
1696     flowtrack->SetSource(AliFlowTrack::kFromKink);
1697   }
1698
1699   flowtrack->SetMass(fTrackMass);
1700
1701   return flowtrack;
1702 }
1703
1704 //-----------------------------------------------------------------------
1705 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackGeneric(TObjArray* trackCollection, Int_t trackIndex) const
1706 {
1707   //fill a flow track from tracklet,vzero,pmd,...
1708   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1709   if (flowtrack)
1710   {
1711     flowtrack->Clear();
1712   }
1713   else
1714   {
1715     flowtrack = new AliFlowTrack();
1716     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1717   }
1718
1719   if (FillFlowTrackGeneric(flowtrack)) return flowtrack;
1720   else
1721   {
1722     trackCollection->RemoveAt(trackIndex);
1723     delete flowtrack;
1724     return NULL;
1725   }
1726 }
1727
1728 //-----------------------------------------------------------------------
1729 Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
1730 {
1731   //fill a flow track from tracklet,vzero,pmd,...
1732   TParticle *tmpTParticle=NULL;
1733   AliMCParticle* tmpAliMCParticle=NULL;
1734   switch (fParamMix)
1735   {
1736     case kPure:
1737       flowtrack->SetPhi(fTrackPhi);
1738       flowtrack->SetEta(fTrackEta);
1739       flowtrack->SetWeight(fTrackWeight);
1740       flowtrack->SetPt(fTrackPt);
1741       flowtrack->SetMass(fTrackMass);
1742       break;
1743     case kTrackWithMCkine:
1744       if (!fMCparticle) return kFALSE;
1745       flowtrack->SetPhi( fMCparticle->Phi() );
1746       flowtrack->SetEta( fMCparticle->Eta() );
1747       flowtrack->SetPt( fMCparticle->Pt() );
1748       flowtrack->SetMass(fTrackMass);
1749       break;
1750     case kTrackWithMCpt:
1751       if (!fMCparticle) return kFALSE;
1752       flowtrack->SetPhi(fTrackPhi);
1753       flowtrack->SetEta(fTrackEta);
1754       flowtrack->SetWeight(fTrackWeight);
1755       flowtrack->SetPt(fMCparticle->Pt());
1756       flowtrack->SetMass(fTrackMass);
1757       break;
1758     case kTrackWithPtFromFirstMother:
1759       if (!fMCparticle) return kFALSE;
1760       flowtrack->SetPhi(fTrackPhi);
1761       flowtrack->SetEta(fTrackEta);
1762       flowtrack->SetWeight(fTrackWeight);
1763       tmpTParticle = fMCparticle->Particle();
1764       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1765       flowtrack->SetPt(tmpAliMCParticle->Pt());
1766       flowtrack->SetMass(fTrackMass);
1767       break;
1768     default:
1769       flowtrack->SetPhi(fTrackPhi);
1770       flowtrack->SetEta(fTrackEta);
1771       flowtrack->SetWeight(fTrackWeight);
1772       flowtrack->SetPt(fTrackPt);
1773       flowtrack->SetMass(fTrackMass);
1774       break;
1775   }
1776   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1777   return kTRUE;
1778 }
1779
1780 //-----------------------------------------------------------------------
1781 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackVParticle(TObjArray* trackCollection, Int_t trackIndex) const
1782 {
1783   //fill flow track from AliVParticle (ESD,AOD,MC)
1784
1785   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1786   if (flowtrack)
1787   {
1788     flowtrack->Clear();
1789   }
1790   else
1791   {
1792     flowtrack = new AliFlowTrack();
1793     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1794   }
1795
1796   if (FillFlowTrackVParticle(flowtrack)) return flowtrack;
1797   else
1798   {
1799     trackCollection->RemoveAt(trackIndex);
1800     delete flowtrack;
1801     return NULL;
1802   }
1803 }
1804
1805 //-----------------------------------------------------------------------
1806 Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
1807 {
1808   //fill flow track from AliVParticle (ESD,AOD,MC)
1809   if (!fTrack) return kFALSE;
1810   TParticle *tmpTParticle=NULL;
1811   AliMCParticle* tmpAliMCParticle=NULL;
1812   AliExternalTrackParam* externalParams=NULL;
1813   AliESDtrack* esdtrack=NULL;
1814   switch(fParamMix)
1815   {
1816     case kPure:
1817       flowtrack->Set(fTrack);
1818       break;
1819     case kTrackWithMCkine:
1820       flowtrack->Set(fMCparticle);
1821       break;
1822     case kTrackWithMCPID:
1823       flowtrack->Set(fTrack);
1824       //flowtrack->setPID(...) from mc, when implemented
1825       break;
1826     case kTrackWithMCpt:
1827       if (!fMCparticle) return kFALSE;
1828       flowtrack->Set(fTrack);
1829       flowtrack->SetPt(fMCparticle->Pt());
1830       break;
1831     case kTrackWithPtFromFirstMother:
1832       if (!fMCparticle) return kFALSE;
1833       flowtrack->Set(fTrack);
1834       tmpTParticle = fMCparticle->Particle();
1835       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1836       flowtrack->SetPt(tmpAliMCParticle->Pt());
1837       break;
1838     case kTrackWithTPCInnerParams:
1839       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1840       if (!esdtrack) return kFALSE;
1841       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1842       if (!externalParams) return kFALSE;
1843       flowtrack->Set(externalParams);
1844       break;
1845     default:
1846       flowtrack->Set(fTrack);
1847       break;
1848   }
1849   if (fParamType==kMC)
1850   {
1851     flowtrack->SetSource(AliFlowTrack::kFromMC);
1852     flowtrack->SetID(fTrackLabel);
1853   }
1854   else if (dynamic_cast<AliESDtrack*>(fTrack))
1855   {
1856     flowtrack->SetSource(AliFlowTrack::kFromESD);
1857     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1858   }
1859   else if (dynamic_cast<AliESDMuonTrack*>(fTrack))                                  // XZhang 20120604
1860   {                                                                                 // XZhang 20120604
1861     flowtrack->SetSource(AliFlowTrack::kFromMUON);                                  // XZhang 20120604
1862     flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID());  // XZhang 20120604
1863   }                                                                                 // XZhang 20120604
1864   else if (dynamic_cast<AliAODTrack*>(fTrack))
1865   {
1866     if (fParamType==kMUON)                            // XZhang 20120604
1867       flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1868     else                                              // XZhang 20120604
1869       flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1870     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1871   }
1872   else if (dynamic_cast<AliMCParticle*>(fTrack))
1873   {
1874     flowtrack->SetSource(AliFlowTrack::kFromMC);
1875     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1876   }
1877   flowtrack->SetMass(fTrackMass);
1878   return kTRUE;
1879 }
1880
1881 //-----------------------------------------------------------------------
1882 AliFlowTrack* AliFlowTrackCuts::FillFlowTrack(TObjArray* trackCollection, Int_t trackIndex) const
1883 {
1884   //fill a flow track constructed from whatever we applied cuts on
1885   //return true on success
1886   switch (fParamType)
1887   {
1888     case kSPDtracklet:
1889       return FillFlowTrackGeneric(trackCollection, trackIndex);
1890     case kPMD:
1891       return FillFlowTrackGeneric(trackCollection, trackIndex);
1892     case kVZERO:
1893       return FillFlowTrackGeneric(trackCollection, trackIndex);
1894     case kKink:
1895       return FillFlowTrackKink(trackCollection, trackIndex);
1896     //case kV0:
1897     //  return FillFlowTrackV0(trackCollection, trackIndex);
1898     default:
1899       return FillFlowTrackVParticle(trackCollection, trackIndex);
1900   }
1901 }
1902
1903 //-----------------------------------------------------------------------
1904 Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
1905 {
1906   //fill a flow track constructed from whatever we applied cuts on
1907   //return true on success
1908   switch (fParamType)
1909   {
1910     case kSPDtracklet:
1911       return FillFlowTrackGeneric(track);
1912     case kPMD:
1913       return FillFlowTrackGeneric(track);
1914     case kVZERO:
1915       return FillFlowTrackGeneric(track);
1916     default:
1917       return FillFlowTrackVParticle(track);
1918   }
1919 }
1920
1921 ////-----------------------------------------------------------------------
1922 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
1923 //{
1924 //  //make a flow track from tracklet
1925 //  AliFlowTrack* flowtrack=NULL;
1926 //  TParticle *tmpTParticle=NULL;
1927 //  AliMCParticle* tmpAliMCParticle=NULL;
1928 //  switch (fParamMix)
1929 //  {
1930 //    case kPure:
1931 //      flowtrack = new AliFlowTrack();
1932 //      flowtrack->SetPhi(fTrackPhi);
1933 //      flowtrack->SetEta(fTrackEta);
1934 //      flowtrack->SetWeight(fTrackWeight);
1935 //      break;
1936 //    case kTrackWithMCkine:
1937 //      if (!fMCparticle) return NULL;
1938 //      flowtrack = new AliFlowTrack();
1939 //      flowtrack->SetPhi( fMCparticle->Phi() );
1940 //      flowtrack->SetEta( fMCparticle->Eta() );
1941 //      flowtrack->SetPt( fMCparticle->Pt() );
1942 //      break;
1943 //    case kTrackWithMCpt:
1944 //      if (!fMCparticle) return NULL;
1945 //      flowtrack = new AliFlowTrack();
1946 //      flowtrack->SetPhi(fTrackPhi);
1947 //      flowtrack->SetEta(fTrackEta);
1948 //      flowtrack->SetWeight(fTrackWeight);
1949 //      flowtrack->SetPt(fMCparticle->Pt());
1950 //      break;
1951 //    case kTrackWithPtFromFirstMother:
1952 //      if (!fMCparticle) return NULL;
1953 //      flowtrack = new AliFlowTrack();
1954 //      flowtrack->SetPhi(fTrackPhi);
1955 //      flowtrack->SetEta(fTrackEta);
1956 //      flowtrack->SetWeight(fTrackWeight);
1957 //      tmpTParticle = fMCparticle->Particle();
1958 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1959 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
1960 //      break;
1961 //    default:
1962 //      flowtrack = new AliFlowTrack();
1963 //      flowtrack->SetPhi(fTrackPhi);
1964 //      flowtrack->SetEta(fTrackEta);
1965 //      flowtrack->SetWeight(fTrackWeight);
1966 //      break;
1967 //  }
1968 //  flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1969 //  return flowtrack;
1970 //}
1971 //
1972 ////-----------------------------------------------------------------------
1973 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
1974 //{
1975 //  //make flow track from AliVParticle (ESD,AOD,MC)
1976 //  if (!fTrack) return NULL;
1977 //  AliFlowTrack* flowtrack=NULL;
1978 //  TParticle *tmpTParticle=NULL;
1979 //  AliMCParticle* tmpAliMCParticle=NULL;
1980 //  AliExternalTrackParam* externalParams=NULL;
1981 //  AliESDtrack* esdtrack=NULL;
1982 //  switch(fParamMix)
1983 //  {
1984 //    case kPure:
1985 //      flowtrack = new AliFlowTrack(fTrack);
1986 //      break;
1987 //    case kTrackWithMCkine:
1988 //      flowtrack = new AliFlowTrack(fMCparticle);
1989 //      break;
1990 //    case kTrackWithMCPID:
1991 //      flowtrack = new AliFlowTrack(fTrack);
1992 //      //flowtrack->setPID(...) from mc, when implemented
1993 //      break;
1994 //    case kTrackWithMCpt:
1995 //      if (!fMCparticle) return NULL;
1996 //      flowtrack = new AliFlowTrack(fTrack);
1997 //      flowtrack->SetPt(fMCparticle->Pt());
1998 //      break;
1999 //    case kTrackWithPtFromFirstMother:
2000 //      if (!fMCparticle) return NULL;
2001 //      flowtrack = new AliFlowTrack(fTrack);
2002 //      tmpTParticle = fMCparticle->Particle();
2003 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2004 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2005 //      break;
2006 //    case kTrackWithTPCInnerParams:
2007 //      esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2008 //      if (!esdtrack) return NULL;
2009 //      externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2010 //      if (!externalParams) return NULL;
2011 //      flowtrack = new AliFlowTrack(externalParams);
2012 //      break;
2013 //    default:
2014 //      flowtrack = new AliFlowTrack(fTrack);
2015 //      break;
2016 //  }
2017 //  if (fParamType==kMC)
2018 //  {
2019 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
2020 //    flowtrack->SetID(fTrackLabel);
2021 //  }
2022 //  else if (dynamic_cast<AliESDtrack*>(fTrack))
2023 //  {
2024 //    flowtrack->SetSource(AliFlowTrack::kFromESD);
2025 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2026 //  }
2027 //  else if (dynamic_cast<AliAODTrack*>(fTrack))
2028 //  {
2029 //    flowtrack->SetSource(AliFlowTrack::kFromAOD);
2030 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2031 //  }
2032 //  else if (dynamic_cast<AliMCParticle*>(fTrack))
2033 //  {
2034 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
2035 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2036 //  }
2037 //  return flowtrack;
2038 //}
2039 //
2040 ////-----------------------------------------------------------------------
2041 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
2042 //{
2043 //  //make a flow track from PMD track
2044 //  AliFlowTrack* flowtrack=NULL;
2045 //  TParticle *tmpTParticle=NULL;
2046 //  AliMCParticle* tmpAliMCParticle=NULL;
2047 //  switch (fParamMix)
2048 //  {
2049 //    case kPure:
2050 //      flowtrack = new AliFlowTrack();
2051 //      flowtrack->SetPhi(fTrackPhi);
2052 //      flowtrack->SetEta(fTrackEta);
2053 //      flowtrack->SetWeight(fTrackWeight);
2054 //      break;
2055 //    case kTrackWithMCkine:
2056 //      if (!fMCparticle) return NULL;
2057 //      flowtrack = new AliFlowTrack();
2058 //      flowtrack->SetPhi( fMCparticle->Phi() );
2059 //      flowtrack->SetEta( fMCparticle->Eta() );
2060 //      flowtrack->SetWeight(fTrackWeight);
2061 //      flowtrack->SetPt( fMCparticle->Pt() );
2062 //      break;
2063 //    case kTrackWithMCpt:
2064 //      if (!fMCparticle) return NULL;
2065 //      flowtrack = new AliFlowTrack();
2066 //      flowtrack->SetPhi(fTrackPhi);
2067 //      flowtrack->SetEta(fTrackEta);
2068 //      flowtrack->SetWeight(fTrackWeight);
2069 //      flowtrack->SetPt(fMCparticle->Pt());
2070 //      break;
2071 //    case kTrackWithPtFromFirstMother:
2072 //      if (!fMCparticle) return NULL;
2073 //      flowtrack = new AliFlowTrack();
2074 //      flowtrack->SetPhi(fTrackPhi);
2075 //      flowtrack->SetEta(fTrackEta);
2076 //      flowtrack->SetWeight(fTrackWeight);
2077 //      tmpTParticle = fMCparticle->Particle();
2078 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2079 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2080 //      break;
2081 //    default:
2082 //      flowtrack = new AliFlowTrack();
2083 //      flowtrack->SetPhi(fTrackPhi);
2084 //      flowtrack->SetEta(fTrackEta);
2085 //      flowtrack->SetWeight(fTrackWeight);
2086 //      break;
2087 //  }
2088 //
2089 //  flowtrack->SetSource(AliFlowTrack::kFromPMD);
2090 //  return flowtrack;
2091 //}
2092 //
2093 ////-----------------------------------------------------------------------
2094 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVZERO() const
2095 //{
2096 //  //make a flow track from VZERO
2097 //  AliFlowTrack* flowtrack=NULL;
2098 //  TParticle *tmpTParticle=NULL;
2099 //  AliMCParticle* tmpAliMCParticle=NULL;
2100 //  switch (fParamMix)
2101 //  {
2102 //    case kPure:
2103 //      flowtrack = new AliFlowTrack();
2104 //      flowtrack->SetPhi(fTrackPhi);
2105 //      flowtrack->SetEta(fTrackEta);
2106 //      flowtrack->SetWeight(fTrackWeight);
2107 //      break;
2108 //    case kTrackWithMCkine:
2109 //      if (!fMCparticle) return NULL;
2110 //      flowtrack = new AliFlowTrack();
2111 //      flowtrack->SetPhi( fMCparticle->Phi() );
2112 //      flowtrack->SetEta( fMCparticle->Eta() );
2113 //      flowtrack->SetWeight(fTrackWeight);
2114 //      flowtrack->SetPt( fMCparticle->Pt() );
2115 //      break;
2116 //    case kTrackWithMCpt:
2117 //      if (!fMCparticle) return NULL;
2118 //      flowtrack = new AliFlowTrack();
2119 //      flowtrack->SetPhi(fTrackPhi);
2120 //      flowtrack->SetEta(fTrackEta);
2121 //      flowtrack->SetWeight(fTrackWeight);
2122 //      flowtrack->SetPt(fMCparticle->Pt());
2123 //      break;
2124 //    case kTrackWithPtFromFirstMother:
2125 //      if (!fMCparticle) return NULL;
2126 //      flowtrack = new AliFlowTrack();
2127 //      flowtrack->SetPhi(fTrackPhi);
2128 //      flowtrack->SetEta(fTrackEta);
2129 //      flowtrack->SetWeight(fTrackWeight);
2130 //      tmpTParticle = fMCparticle->Particle();
2131 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2132 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2133 //      break;
2134 //    default:
2135 //      flowtrack = new AliFlowTrack();
2136 //      flowtrack->SetPhi(fTrackPhi);
2137 //      flowtrack->SetEta(fTrackEta);
2138 //      flowtrack->SetWeight(fTrackWeight);
2139 //      break;
2140 //  }
2141 //
2142 //  flowtrack->SetSource(AliFlowTrack::kFromVZERO);
2143 //  return flowtrack;
2144 //}
2145 //
2146 ////-----------------------------------------------------------------------
2147 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
2148 //{
2149 //  //get a flow track constructed from whatever we applied cuts on
2150 //  //caller is resposible for deletion
2151 //  //if construction fails return NULL
2152 //  //TODO: for tracklets, PMD and VZERO we probably need just one method,
2153 //  //something like MakeFlowTrackGeneric(), wait with this until
2154 //  //requirements quirks are known.
2155 //  switch (fParamType)
2156 //  {
2157 //    case kSPDtracklet:
2158 //      return MakeFlowTrackSPDtracklet();
2159 //    case kPMD:
2160 //      return MakeFlowTrackPMDtrack();
2161 //    case kVZERO:
2162 //      return MakeFlowTrackVZERO();
2163 //    default:
2164 //      return MakeFlowTrackVParticle();
2165 //  }
2166 //}
2167
2168 //-----------------------------------------------------------------------
2169 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
2170 {
2171   //check if current particle is a physical primary
2172   if (!fMCevent) return kFALSE;
2173   if (fTrackLabel<0) return kFALSE;
2174   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
2175 }
2176
2177 //-----------------------------------------------------------------------
2178 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
2179 {
2180   //check if current particle is a physical primary
2181   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
2182   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
2183   if (!track) return kFALSE;
2184   TParticle* particle = track->Particle();
2185   Bool_t transported = particle->TestBit(kTransportBit);
2186   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
2187         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
2188   return (physprim && (transported || !requiretransported));
2189 }
2190
2191 //-----------------------------------------------------------------------
2192 void AliFlowTrackCuts::DefineHistograms()
2193 {
2194   //define qa histograms
2195   if (fQA) return;
2196
2197   const Int_t kNbinsP=200;
2198   Double_t binsP[kNbinsP+1];
2199   binsP[0]=0.0;
2200   for(int i=1; i<kNbinsP+1; i++)
2201   {
2202     //if(binsP[i-1]+0.05<1.01)
2203     //  binsP[i]=binsP[i-1]+0.05;
2204     //else
2205     binsP[i]=binsP[i-1]+0.05;
2206   }
2207
2208   const Int_t nBinsDCA=1000;
2209   Double_t binsDCA[nBinsDCA+1];
2210   for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
2211   //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
2212
2213   Bool_t adddirstatus = TH1::AddDirectoryStatus();
2214   TH1::AddDirectory(kFALSE);
2215   fQA=new TList(); fQA->SetOwner();
2216   fQA->SetName(Form("%s QA",GetName()));
2217   TList* before = new TList(); before->SetOwner();
2218   before->SetName("before");
2219   TList* after = new TList(); after->SetOwner();
2220   after->SetName("after");
2221   fQA->Add(before);
2222   fQA->Add(after);
2223   before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2224   after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2225   before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2226   after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2227   before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2228   after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2229   //primary
2230   TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2231   TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2232   TAxis* axis = NULL;
2233   axis = hb->GetYaxis();
2234   axis->SetBinLabel(1,"secondary");
2235   axis->SetBinLabel(2,"primary");
2236   axis = ha->GetYaxis();
2237   axis->SetBinLabel(1,"secondary");
2238   axis->SetBinLabel(2,"primary");
2239   before->Add(hb); //3
2240   after->Add(ha); //3
2241   //production process
2242   hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2243                       -0.5, kMaxMCProcess-0.5);
2244   ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2245                       -0.5, kMaxMCProcess-0.5);
2246   axis = hb->GetYaxis();
2247   for (Int_t i=0; i<kMaxMCProcess; i++)
2248   {
2249     axis->SetBinLabel(i+1,TMCProcessName[i]);
2250   }
2251   axis = ha->GetYaxis();
2252   for (Int_t i=0; i<kMaxMCProcess; i++)
2253   {
2254     axis->SetBinLabel(i+1,TMCProcessName[i]);
2255   }
2256   before->Add(hb); //4
2257   after->Add(ha); //4
2258   //DCA
2259   before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2260   after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2261   before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2262   after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2263   //first mother
2264   hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2265   ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2266   hb->GetYaxis()->SetBinLabel(1, "#gamma");
2267   ha->GetYaxis()->SetBinLabel(1, "#gamma");
2268   hb->GetYaxis()->SetBinLabel(2, "e^{+}");
2269   ha->GetYaxis()->SetBinLabel(2, "e^{+}");
2270   hb->GetYaxis()->SetBinLabel(3, "e^{-}");
2271   ha->GetYaxis()->SetBinLabel(3, "e^{-}");
2272   hb->GetYaxis()->SetBinLabel(4, "#nu");
2273   ha->GetYaxis()->SetBinLabel(4, "#nu");
2274   hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2275   ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2276   hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2277   ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2278   hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2279   ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2280   hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2281   ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2282   hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2283   ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2284   hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2285   ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2286   hb->GetYaxis()->SetBinLabel(11, "K^{+}");
2287   ha->GetYaxis()->SetBinLabel(11, "K^{+}");
2288   hb->GetYaxis()->SetBinLabel(12, "K^{-}");
2289   ha->GetYaxis()->SetBinLabel(12, "K^{-}");
2290   hb->GetYaxis()->SetBinLabel( 13, "n");
2291   ha->GetYaxis()->SetBinLabel( 13, "n");
2292   hb->GetYaxis()->SetBinLabel( 14, "p");
2293   ha->GetYaxis()->SetBinLabel( 14, "p");
2294   hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
2295   ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
2296   hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2297   ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2298   hb->GetYaxis()->SetBinLabel(17, "#eta");
2299   ha->GetYaxis()->SetBinLabel(17, "#eta");
2300   hb->GetYaxis()->SetBinLabel(18, "#Lambda");
2301   ha->GetYaxis()->SetBinLabel(18, "#Lambda");
2302   hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2303   ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2304   hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2305   ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2306   hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2307   ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2308   hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2309   ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2310   hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2311   ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2312   hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2313   ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2314   hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
2315   ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
2316   hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2317   ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2318   hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2319   ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2320   hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2321   ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2322   hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2323   ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2324   hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2325   ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2326   hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2327   ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2328   hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2329   ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2330   hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2331   ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2332   hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2333   ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2334   hb->GetYaxis()->SetBinLabel(35, "D^{+}");
2335   ha->GetYaxis()->SetBinLabel(35, "D^{+}");
2336   hb->GetYaxis()->SetBinLabel(36, "D^{-}");
2337   ha->GetYaxis()->SetBinLabel(36, "D^{-}");
2338   hb->GetYaxis()->SetBinLabel(37, "D^{0}");
2339   ha->GetYaxis()->SetBinLabel(37, "D^{0}");
2340   hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2341   ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2342   hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2343   ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2344   hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2345   ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2346   hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2347   ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2348   hb->GetYaxis()->SetBinLabel(42, "W^{+}");
2349   ha->GetYaxis()->SetBinLabel(42, "W^{+}");
2350   hb->GetYaxis()->SetBinLabel(43, "W^{-}");
2351   ha->GetYaxis()->SetBinLabel(43, "W^{-}");
2352   hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
2353   ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
2354   before->Add(hb);//7
2355   after->Add(ha);//7
2356
2357   before->Add(new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2358   after->Add( new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2359
2360   before->Add(new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
2361   after->Add( new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
2362
2363   before->Add(new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
2364   after->Add( new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
2365
2366   before->Add(new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
2367   after->Add( new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
2368
2369   before->Add(new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
2370   after->Add( new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
2371
2372   //kink stuff
2373   before->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2374   after->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2375
2376   before->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2377   after->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2378
2379   before->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2380   after->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2381
2382   before->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2383   after->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2384
2385   before->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2386   after->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2387   TH1::AddDirectory(adddirstatus);
2388 }
2389
2390 //-----------------------------------------------------------------------
2391 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
2392 {
2393   //get the number of tracks in the input event according source
2394   //selection (ESD tracks, tracklets, MC particles etc.)
2395   AliESDEvent* esd=NULL;
2396   AliAODEvent* aod=NULL;  // XZhang 20120615
2397   switch (fParamType)
2398   {
2399     case kSPDtracklet:
2400       if (!fEvent) return 0;                                           // XZhang 20120615
2401       esd = dynamic_cast<AliESDEvent*>(fEvent);
2402       aod = dynamic_cast<AliAODEvent*>(fEvent);                        // XZhang 20120615
2403 //    if (!esd) return 0;                                              // XZhang 20120615
2404 //    return esd->GetMultiplicity()->GetNumberOfTracklets();           // XZhang 20120615
2405       if (esd) return esd->GetMultiplicity()->GetNumberOfTracklets();  // XZhang 20120615
2406       if (aod) return aod->GetTracklets()->GetNumberOfTracklets();     // XZhang 20120615
2407     case kMC:
2408       if (!fMCevent) return 0;
2409       return fMCevent->GetNumberOfTracks();
2410     case kPMD:
2411       esd = dynamic_cast<AliESDEvent*>(fEvent);
2412       if (!esd) return 0;
2413       return esd->GetNumberOfPmdTracks();
2414     case kVZERO:
2415       return fgkNumberOfVZEROtracks;
2416     case kMUON:                                      // XZhang 20120604
2417       if (!fEvent) return 0;                         // XZhang 20120604
2418       esd = dynamic_cast<AliESDEvent*>(fEvent);      // XZhang 20120604
2419       if (esd) return esd->GetNumberOfMuonTracks();  // XZhang 20120604
2420       return fEvent->GetNumberOfTracks();  // if AOD // XZhang 20120604
2421     case kKink:
2422       esd = dynamic_cast<AliESDEvent*>(fEvent);
2423       if (!esd) return 0;
2424       return esd->GetNumberOfKinks();
2425     case kV0:
2426       esd = dynamic_cast<AliESDEvent*>(fEvent);
2427       if (!esd) return 0;
2428       return esd->GetNumberOfV0s();
2429     default:
2430       if (!fEvent) return 0;
2431       return fEvent->GetNumberOfTracks();
2432   }
2433   return 0;
2434 }
2435
2436 //-----------------------------------------------------------------------
2437 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
2438 {
2439   //get the input object according the data source selection:
2440   //(esd tracks, traclets, mc particles,etc...)
2441   AliESDEvent* esd=NULL;
2442   AliAODEvent* aod=NULL;  // XZhang 20120615
2443   switch (fParamType)
2444   {
2445     case kSPDtracklet:
2446       if (!fEvent) return NULL;                                              // XZhang 20120615
2447       esd = dynamic_cast<AliESDEvent*>(fEvent);
2448       aod = dynamic_cast<AliAODEvent*>(fEvent);                              // XZhang 20120615
2449 //    if (!esd) return NULL;                                                 // XZhang 20120615
2450 //    return const_cast<AliMultiplicity*>(esd->GetMultiplicity());           // XZhang 20120615
2451       if (esd) return const_cast<AliMultiplicity*>(esd->GetMultiplicity());  // XZhang 20120615
2452       if (aod) return const_cast<AliAODTracklets*>(aod->GetTracklets());     // XZhang 20120615
2453     case kMC:
2454       if (!fMCevent) return NULL;
2455       return fMCevent->GetTrack(i);
2456     case kPMD:
2457       esd = dynamic_cast<AliESDEvent*>(fEvent);
2458       if (!esd) return NULL;
2459       return esd->GetPmdTrack(i);
2460     case kVZERO:
2461       esd = dynamic_cast<AliESDEvent*>(fEvent);
2462       if (!esd) //contributed by G.Ortona
2463       {
2464         aod = dynamic_cast<AliAODEvent*>(fEvent);
2465         if(!aod)return NULL;
2466         return aod->GetVZEROData();
2467       }
2468       return esd->GetVZEROData();
2469     case kMUON:                                  // XZhang 20120604
2470       if (!fEvent) return NULL;                  // XZhang 20120604
2471       esd = dynamic_cast<AliESDEvent*>(fEvent);  // XZhang 20120604
2472       if (esd) return esd->GetMuonTrack(i);      // XZhang 20120604
2473       return fEvent->GetTrack(i);  // if AOD     // XZhang 20120604
2474     case kKink:
2475       esd = dynamic_cast<AliESDEvent*>(fEvent);
2476       if (!esd) return NULL;
2477       return esd->GetKink(i);
2478     case kV0:
2479       esd = dynamic_cast<AliESDEvent*>(fEvent);
2480       if (!esd) return NULL;
2481       return esd->GetV0(i);
2482     default:
2483       if (!fEvent) return NULL;
2484       return fEvent->GetTrack(i);
2485   }
2486 }
2487
2488 //-----------------------------------------------------------------------
2489 void AliFlowTrackCuts::Clear(Option_t*)
2490 {
2491   //clean up
2492   fMCevent=NULL;
2493   fEvent=NULL;
2494   ClearTrack();
2495 }
2496
2497 //-----------------------------------------------------------------------
2498 void AliFlowTrackCuts::ClearTrack(Option_t*)
2499 {
2500   //clean up last track
2501   fKink=NULL;
2502   fV0=NULL;
2503   fTrack=NULL;
2504   fMCparticle=NULL;
2505   fTrackLabel=-997;
2506   fTrackWeight=1.0;
2507   fTrackEta=0.0;
2508   fTrackPhi=0.0;
2509   fTrackPt=0.0;
2510   fPOItype=1;
2511   fTrackMass=0.;
2512 }
2513 //-----------------------------------------------------------------------
2514 Bool_t AliFlowTrackCuts::PassesAODpidCut(const AliAODTrack* track )
2515 {
2516      if(!track->GetAODEvent()->GetTOFHeader()){
2517           AliAODPid *pidObj = track->GetDetPid();
2518           if (!pidObj) fESDpid.GetTOFResponse().SetTimeResolution(84.);
2519           else{
2520             Double_t sigmaTOFPidInAOD[10];
2521             pidObj->GetTOFpidResolution(sigmaTOFPidInAOD);
2522             if(sigmaTOFPidInAOD[0] > 84.){
2523               fESDpid.GetTOFResponse().SetTimeResolution(sigmaTOFPidInAOD[0]); // use the electron TOF PID sigma as time resolution (including the T0 used)
2524           }
2525         }
2526      }
2527
2528  //check if passes the selected pid cut for ESDs
2529   Bool_t pass = kTRUE;
2530   switch (fPIDsource)
2531   {
2532    case kTOFbeta:
2533       if (!PassesTOFbetaCut(track)) pass=kFALSE;
2534       break;
2535   case kTOFbayesian:
2536       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2537       break;
2538   case kTPCbayesian:
2539       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2540       break;
2541   default:
2542     return kTRUE;
2543     break;
2544  }
2545   return pass;
2546
2547 }
2548 //-----------------------------------------------------------------------
2549 Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
2550 {
2551   //check if passes the selected pid cut for ESDs
2552   Bool_t pass = kTRUE;
2553   switch (fPIDsource)
2554   {
2555     case kTPCpid:
2556       if (!PassesTPCpidCut(track)) pass=kFALSE;
2557       break;
2558     case kTPCdedx:
2559       if (!PassesTPCdedxCut(track)) pass=kFALSE;
2560       break;
2561     case kTOFpid:
2562       if (!PassesTOFpidCut(track)) pass=kFALSE;
2563       break;
2564     case kTOFbeta:
2565       if (!PassesTOFbetaCut(track)) pass=kFALSE;
2566       break;
2567     case kTOFbetaSimple:
2568       if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
2569       break;
2570     case kTPCbayesian:
2571       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2572       break;
2573       // part added by F. Noferini
2574     case kTOFbayesian:
2575       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2576       break;
2577       // end part added by F. Noferini
2578
2579       //part added by Natasha
2580     case kTPCNuclei:
2581       if (!PassesNucleiSelection(track)) pass=kFALSE;
2582       break;
2583       //end part added by Natasha
2584
2585     default:
2586       printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
2587       pass=kFALSE;
2588       break;
2589   }
2590   return pass;
2591 }
2592
2593 //-----------------------------------------------------------------------
2594 Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
2595 {
2596   //check if passes PID cut using timing in TOF
2597   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2598                      (track->GetTOFsignal() > 12000) &&
2599                      (track->GetTOFsignal() < 100000) &&
2600                      (track->GetIntegratedLength() > 365);
2601
2602   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2603
2604   Bool_t statusMatchingHard = TPCTOFagree(track);
2605   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2606        return kFALSE;
2607
2608   if (!goodtrack) return kFALSE;
2609
2610   const Float_t c = 2.99792457999999984e-02;
2611   Float_t p = track->GetP();
2612   Float_t l = track->GetIntegratedLength();
2613   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2614   Float_t timeTOF = track->GetTOFsignal()- trackT0;
2615   Float_t beta = l/timeTOF/c;
2616   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2617   track->GetIntegratedTimes(integratedTimes);
2618   Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
2619   Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
2620   for (Int_t i=0;i<5;i++)
2621   {
2622     betaHypothesis[i] = l/integratedTimes[i]/c;
2623     s[i] = beta-betaHypothesis[i];
2624   }
2625
2626   switch (fParticleID)
2627   {
2628     case AliPID::kPion:
2629       return ( (s[2]<0.015) && (s[2]>-0.015) &&
2630                (s[3]>0.025) &&
2631                (s[4]>0.03) );
2632     case AliPID::kKaon:
2633       return ( (s[3]<0.015) && (s[3]>-0.015) &&
2634                (s[2]<-0.03) &&
2635                (s[4]>0.03) );
2636     case AliPID::kProton:
2637       return ( (s[4]<0.015) && (s[4]>-0.015) &&
2638                (s[3]<-0.025) &&
2639                (s[2]<-0.025) );
2640     default:
2641       return kFALSE;
2642   }
2643   return kFALSE;
2644 }
2645
2646 //-----------------------------------------------------------------------
2647 Float_t AliFlowTrackCuts::GetBeta(const AliVTrack* track)
2648 {
2649   //get beta
2650   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2651   track->GetIntegratedTimes(integratedTimes);
2652
2653   const Float_t c = 2.99792457999999984e-02;
2654   Float_t p = track->P();
2655   Float_t l = integratedTimes[0]*c;
2656   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2657   Float_t timeTOF = track->GetTOFsignal()- trackT0;
2658   return l/timeTOF/c;
2659 }
2660 //-----------------------------------------------------------------------
2661 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliAODTrack* track )
2662 {
2663   //check if track passes pid selection with an asymmetric TOF beta cut
2664   if (!fTOFpidCuts)
2665   {
2666     //printf("no TOFpidCuts\n");
2667     return kFALSE;
2668   }
2669
2670   //check if passes PID cut using timing in TOF
2671   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2672                      (track->GetTOFsignal() > 12000) &&
2673                      (track->GetTOFsignal() < 100000);
2674
2675   if (!goodtrack) return kFALSE;
2676
2677   const Float_t c = 2.99792457999999984e-02;
2678   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2679   track->GetIntegratedTimes(integratedTimes);
2680   Float_t l = integratedTimes[0]*c;
2681
2682   goodtrack = goodtrack && (l > 365);
2683
2684   if (!goodtrack) return kFALSE;
2685
2686   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2687
2688   Bool_t statusMatchingHard = TPCTOFagree(track);
2689   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2690        return kFALSE;
2691
2692
2693   Float_t beta = GetBeta(track);
2694
2695   //construct the pid index because it's not AliPID::EParticleType
2696   Int_t pid = 0;
2697   switch (fParticleID)
2698   {
2699     case AliPID::kPion:
2700       pid=2;
2701       break;
2702     case AliPID::kKaon:
2703       pid=3;
2704       break;
2705     case AliPID::kProton:
2706       pid=4;
2707       break;
2708     default:
2709       return kFALSE;
2710   }
2711
2712   //signal to cut on
2713   Float_t p = track->P();
2714   Float_t betahypothesis = l/integratedTimes[pid]/c;
2715   Float_t betadiff = beta-betahypothesis;
2716
2717   Float_t* arr = fTOFpidCuts->GetMatrixArray();
2718   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2719   if (col<0) return kFALSE;
2720   Float_t min = (*fTOFpidCuts)(1,col);
2721   Float_t max = (*fTOFpidCuts)(2,col);
2722
2723   Bool_t pass = (betadiff>min && betadiff<max);
2724
2725   return pass;
2726 }
2727 //-----------------------------------------------------------------------
2728 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
2729 {
2730   //check if track passes pid selection with an asymmetric TOF beta cut
2731   if (!fTOFpidCuts)
2732   {
2733     //printf("no TOFpidCuts\n");
2734     return kFALSE;
2735   }
2736
2737   //check if passes PID cut using timing in TOF
2738   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2739                      (track->GetTOFsignal() > 12000) &&
2740                      (track->GetTOFsignal() < 100000) &&
2741                      (track->GetIntegratedLength() > 365);
2742
2743   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2744
2745   if (!goodtrack) return kFALSE;
2746
2747   Bool_t statusMatchingHard = TPCTOFagree(track);
2748   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2749        return kFALSE;
2750
2751   Float_t beta = GetBeta(track);
2752   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2753   track->GetIntegratedTimes(integratedTimes);
2754
2755   //construct the pid index because it's not AliPID::EParticleType
2756   Int_t pid = 0;
2757   switch (fParticleID)
2758   {
2759     case AliPID::kPion:
2760       pid=2;
2761       break;
2762     case AliPID::kKaon:
2763       pid=3;
2764       break;
2765     case AliPID::kProton:
2766       pid=4;
2767       break;
2768     default:
2769       return kFALSE;
2770   }
2771
2772   //signal to cut on
2773   const Float_t c = 2.99792457999999984e-02;
2774   Float_t l = track->GetIntegratedLength();
2775   Float_t p = track->GetP();
2776   Float_t betahypothesis = l/integratedTimes[pid]/c;
2777   Float_t betadiff = beta-betahypothesis;
2778
2779   Float_t* arr = fTOFpidCuts->GetMatrixArray();
2780   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2781   if (col<0) return kFALSE;
2782   Float_t min = (*fTOFpidCuts)(1,col);
2783   Float_t max = (*fTOFpidCuts)(2,col);
2784
2785   Bool_t pass = (betadiff>min && betadiff<max);
2786
2787   return pass;
2788 }
2789
2790 //-----------------------------------------------------------------------
2791 Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
2792 {
2793   //check if passes PID cut using default TOF pid
2794   Double_t pidTOF[AliPID::kSPECIES];
2795   track->GetTOFpid(pidTOF);
2796   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
2797   return kFALSE;
2798 }
2799
2800 //-----------------------------------------------------------------------
2801 Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
2802 {
2803   //check if passes PID cut using default TPC pid
2804   Double_t pidTPC[AliPID::kSPECIES];
2805   track->GetTPCpid(pidTPC);
2806   Double_t probablity = 0.;
2807   switch (fParticleID)
2808   {
2809     case AliPID::kPion:
2810       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
2811       break;
2812     default:
2813       probablity = pidTPC[fParticleID];
2814   }
2815   if (probablity >= fParticleProbability) return kTRUE;
2816   return kFALSE;
2817 }
2818
2819 //-----------------------------------------------------------------------
2820 Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
2821 {
2822   //get TPC dedx
2823   return track->GetTPCsignal();
2824 }
2825
2826 //-----------------------------------------------------------------------
2827 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
2828 {
2829   //check if passes PID cut using dedx signal in the TPC
2830   if (!fTPCpidCuts)
2831   {
2832     //printf("no TPCpidCuts\n");
2833     return kFALSE;
2834   }
2835
2836   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
2837   if (!tpcparam) return kFALSE;
2838   Double_t p = tpcparam->GetP();
2839   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
2840   Float_t sigTPC = track->GetTPCsignal();
2841   Float_t s = (sigTPC-sigExp)/sigExp;
2842
2843   Float_t* arr = fTPCpidCuts->GetMatrixArray();
2844   Int_t arrSize = fTPCpidCuts->GetNcols();
2845   Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
2846   if (col<0) return kFALSE;
2847   Float_t min = (*fTPCpidCuts)(1,col);
2848   Float_t max = (*fTPCpidCuts)(2,col);
2849
2850   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
2851   return (s>min && s<max);
2852 }
2853
2854 //-----------------------------------------------------------------------
2855 void AliFlowTrackCuts::InitPIDcuts()
2856 {
2857   //init matrices with PID cuts
2858   TMatrixF* t = NULL;
2859   if (!fTPCpidCuts)
2860   {
2861     if (fParticleID==AliPID::kPion)
2862     {
2863       t = new TMatrixF(3,15);
2864       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.0;
2865       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.1;
2866       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.2;
2867       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.2;
2868       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
2869       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
2870       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
2871       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
2872       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
2873       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.4;  (*t)(2,9)  =  0.05;
2874       (*t)(0,10)  = 0.70;  (*t)(1,10)  = -0.4;  (*t)(2,10)  =     0;
2875       (*t)(0,11)  = 0.75;  (*t)(1,11)  = -0.4;  (*t)(2,11)  =     0;
2876       (*t)(0,12)  = 0.80;  (*t)(1,12)  = -0.4;  (*t)(2,12)  = -0.05;
2877       (*t)(0,13)  = 0.85;  (*t)(1,13)  = -0.4;  (*t)(2,13)  = -0.1;
2878       (*t)(0,14)  = 0.90;  (*t)(1,14)  = 0;     (*t)(2,14)  =     0;
2879     }
2880     else
2881     if (fParticleID==AliPID::kKaon)
2882     {
2883       t = new TMatrixF(3,12);
2884       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.2;
2885       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  = 0.2;
2886       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.2;  (*t)(2,2)  = 0.2;
2887       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.2;  (*t)(2,3)  = 0.2;
2888       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.2;
2889       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.2;
2890       (*t)(0,6)  = 0.50;  (*t)(1,6)  =-0.05;  (*t)(2,6)  = 0.2;
2891       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.1;  (*t)(2,7)  = 0.1;
2892       (*t)(0,8)  = 0.60;  (*t)(1,8)  =-0.05;  (*t)(2,8)  = 0.1;
2893       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  = 0.15;
2894       (*t)(0,10)  = 0.70;  (*t)(1,10)  = 0.05;  (*t)(2,10)  = 0.2;
2895       (*t)(0,11)  = 0.75;  (*t)(1,11)  =    0;  (*t)(2,11)  = 0;
2896     }
2897     else
2898     if (fParticleID==AliPID::kProton)
2899     {
2900       t = new TMatrixF(3,9);
2901       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.1;  (*t)(2,0)  =  0.1;
2902       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  =  0.2;
2903       (*t)(0,2)  = 0.80;  (*t)(1,2)  = -0.1;  (*t)(2,2)  =  0.2;
2904       (*t)(0,3)  = 0.85;  (*t)(1,3)  =-0.05;  (*t)(2,3)  =  0.2;
2905       (*t)(0,4)  = 0.90;  (*t)(1,4)  =-0.05;  (*t)(2,4)  = 0.25;
2906       (*t)(0,5)  = 0.95;  (*t)(1,5)  =-0.05;  (*t)(2,5)  = 0.25;
2907       (*t)(0,6)  = 1.00;  (*t)(1,6)  = -0.1;  (*t)(2,6)  = 0.25;
2908       (*t)(0,7)  = 1.10;  (*t)(1,7)  =-0.05;  (*t)(2,7)  =  0.3;
2909       (*t)(0,8) = 1.20;   (*t)(1,8)  =    0;  (*t)(2,8) =    0;
2910     }
2911     delete fTPCpidCuts;
2912     fTPCpidCuts=t;
2913   }
2914   t = NULL;
2915   if (!fTOFpidCuts)
2916   {
2917     if (fParticleID==AliPID::kPion)
2918     {
2919       //TOF pions, 0.9 purity
2920       t = new TMatrixF(3,61);
2921       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2922       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2923       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2924       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2925       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.030;  (*t)(2,4)  =   0.030;
2926       (*t)(0,5)  = 0.250;  (*t)(1,5)  = -0.036;  (*t)(2,5)  =   0.032;
2927       (*t)(0,6)  = 0.300;  (*t)(1,6)  = -0.038;  (*t)(2,6)  =   0.032;
2928       (*t)(0,7)  = 0.350;  (*t)(1,7)  = -0.034;  (*t)(2,7)  =   0.032;
2929       (*t)(0,8)  = 0.400;  (*t)(1,8)  = -0.032;  (*t)(2,8)  =   0.020;
2930       (*t)(0,9)  = 0.450;  (*t)(1,9)  = -0.030;  (*t)(2,9)  =   0.020;
2931       (*t)(0,10)  = 0.500;  (*t)(1,10)  = -0.030;  (*t)(2,10)  =   0.020;
2932       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.030;  (*t)(2,11)  =   0.020;
2933       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.030;  (*t)(2,12)  =   0.020;
2934       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.030;  (*t)(2,13)  =   0.020;
2935       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.030;  (*t)(2,14)  =   0.020;
2936       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.030;  (*t)(2,15)  =   0.020;
2937       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.030;  (*t)(2,16)  =   0.020;
2938       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.030;  (*t)(2,17)  =   0.020;
2939       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.030;  (*t)(2,18)  =   0.020;
2940       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.028;  (*t)(2,19)  =   0.028;
2941       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.028;  (*t)(2,20)  =   0.028;
2942       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.028;  (*t)(2,21)  =   0.028;
2943       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.028;
2944       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.028;
2945       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.020;  (*t)(2,24)  =   0.028;
2946       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.018;  (*t)(2,25)  =   0.028;
2947       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.016;  (*t)(2,26)  =   0.028;
2948       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.014;  (*t)(2,27)  =   0.028;
2949       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.012;  (*t)(2,28)  =   0.026;
2950       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.010;  (*t)(2,29)  =   0.026;
2951       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.008;  (*t)(2,30)  =   0.026;
2952       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.008;  (*t)(2,31)  =   0.024;
2953       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.006;  (*t)(2,32)  =   0.024;
2954       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.004;  (*t)(2,33)  =   0.024;
2955       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.004;  (*t)(2,34)  =   0.024;
2956       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.002;  (*t)(2,35)  =   0.024;
2957       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.002;  (*t)(2,36)  =   0.024;
2958       (*t)(0,37)  = 2.700;  (*t)(1,37)  = 0.000;  (*t)(2,37)  =   0.024;
2959       (*t)(0,38)  = 2.800;  (*t)(1,38)  = 0.000;  (*t)(2,38)  =   0.026;
2960       (*t)(0,39)  = 2.900;  (*t)(1,39)  = 0.000;  (*t)(2,39)  =   0.024;
2961       (*t)(0,40)  = 3.000;  (*t)(1,40)  = 0.002;  (*t)(2,40)  =   0.026;
2962       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.002;  (*t)(2,41)  =   0.026;
2963       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.002;  (*t)(2,42)  =   0.026;
2964       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.002;  (*t)(2,43)  =   0.026;
2965       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.002;  (*t)(2,44)  =   0.026;
2966       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.002;  (*t)(2,45)  =   0.026;
2967       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.002;  (*t)(2,46)  =   0.026;
2968       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.002;  (*t)(2,47)  =   0.026;
2969       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.002;  (*t)(2,48)  =   0.026;
2970       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.004;  (*t)(2,49)  =   0.024;
2971       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.004;  (*t)(2,50)  =   0.026;
2972       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.004;  (*t)(2,51)  =   0.026;
2973       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.004;  (*t)(2,52)  =   0.024;
2974       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.006;  (*t)(2,53)  =   0.024;
2975       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2976       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2977       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2978       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2979       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2980       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2981       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2982     }
2983     else
2984     if (fParticleID==AliPID::kProton)
2985     {
2986       //TOF protons, 0.9 purity
2987       t = new TMatrixF(3,61);
2988       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2989       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2990       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2991       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2992       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.07;  (*t)(2,4)  =   0.07;
2993       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.07;  (*t)(2,5)  =   0.07;
2994       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.07;  (*t)(2,6)  =   0.07;
2995       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.07;  (*t)(2,7)  =   0.07;
2996       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.07;  (*t)(2,8)  =   0.07;
2997       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.07;  (*t)(2,9)  =   0.07;
2998       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.07;  (*t)(2,10)  =   0.07;
2999       (*t)(0,11)  = 0.200;  (*t)(1,11)  = -0.07;  (*t)(2,11)  =   0.07;
3000       (*t)(0,12)  = 0.200;  (*t)(1,12)  = -0.07;  (*t)(2,12)  =   0.07;
3001       (*t)(0,13)  = 0.200;  (*t)(1,13)  = -0.07;  (*t)(2,13)  =   0.07;
3002       (*t)(0,14)  = 0.200;  (*t)(1,14)  = -0.07;  (*t)(2,14)  =   0.07;
3003       (*t)(0,15)  = 0.200;  (*t)(1,15)  = -0.07;  (*t)(2,15)  =   0.07;
3004       (*t)(0,16)  = 0.200;  (*t)(1,16)  = -0.07;  (*t)(2,16)  =   0.07;
3005       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.070;  (*t)(2,17)  =   0.070;
3006       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.072;  (*t)(2,18)  =   0.072;
3007       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.072;  (*t)(2,19)  =   0.072;
3008       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.074;  (*t)(2,20)  =   0.074;
3009       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.032;  (*t)(2,21)  =   0.032;
3010       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.026;
3011       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.026;  (*t)(2,23)  =   0.026;
3012       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.024;  (*t)(2,24)  =   0.024;
3013       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.024;  (*t)(2,25)  =   0.024;
3014       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.026;  (*t)(2,26)  =   0.026;
3015       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.026;  (*t)(2,27)  =   0.026;
3016       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.026;  (*t)(2,28)  =   0.026;
3017       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.026;  (*t)(2,29)  =   0.026;
3018       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.026;  (*t)(2,30)  =   0.026;
3019       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.026;
3020       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.026;  (*t)(2,32)  =   0.024;
3021       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.028;  (*t)(2,33)  =   0.022;
3022       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.028;  (*t)(2,34)  =   0.020;
3023       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.028;  (*t)(2,35)  =   0.018;
3024       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.028;  (*t)(2,36)  =   0.016;
3025       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.028;  (*t)(2,37)  =   0.016;
3026       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.030;  (*t)(2,38)  =   0.014;
3027       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.030;  (*t)(2,39)  =   0.012;
3028       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.030;  (*t)(2,40)  =   0.012;
3029       (*t)(0,41)  = 3.100;  (*t)(1,41)  = -0.030;  (*t)(2,41)  =   0.010;
3030       (*t)(0,42)  = 3.200;  (*t)(1,42)  = -0.030;  (*t)(2,42)  =   0.010;
3031       (*t)(0,43)  = 3.300;  (*t)(1,43)  = -0.030;  (*t)(2,43)  =   0.010;
3032       (*t)(0,44)  = 3.400;  (*t)(1,44)  = -0.030;  (*t)(2,44)  =   0.008;
3033       (*t)(0,45)  = 3.500;  (*t)(1,45)  = -0.030;  (*t)(2,45)  =   0.008;
3034       (*t)(0,46)  = 3.600;  (*t)(1,46)  = -0.030;  (*t)(2,46)  =   0.008;
3035       (*t)(0,47)  = 3.700;  (*t)(1,47)  = -0.030;  (*t)(2,47)  =   0.006;
3036       (*t)(0,48)  = 3.800;  (*t)(1,48)  = -0.030;  (*t)(2,48)  =   0.006;
3037       (*t)(0,49)  = 3.900;  (*t)(1,49)  = -0.030;  (*t)(2,49)  =   0.006;
3038       (*t)(0,50)  = 4.000;  (*t)(1,50)  = -0.028;  (*t)(2,50)  =   0.004;
3039       (*t)(0,51)  = 4.100;  (*t)(1,51)  = -0.030;  (*t)(2,51)  =   0.004;
3040       (*t)(0,52)  = 4.200;  (*t)(1,52)  = -0.030;  (*t)(2,52)  =   0.004;
3041       (*t)(0,53)  = 4.300;  (*t)(1,53)  = -0.028;  (*t)(2,53)  =   0.002;
3042       (*t)(0,54)  = 4.400;  (*t)(1,54)  = -0.030;  (*t)(2,54)  =   0.002;
3043       (*t)(0,55)  = 4.500;  (*t)(1,55)  = -0.028;  (*t)(2,55)  =   0.002;
3044       (*t)(0,56)  = 4.600;  (*t)(1,56)  = -0.028;  (*t)(2,56)  =   0.002;
3045       (*t)(0,57)  = 4.700;  (*t)(1,57)  = -0.028;  (*t)(2,57)  =   0.000;
3046       (*t)(0,58)  = 4.800;  (*t)(1,58)  = -0.028;  (*t)(2,58)  =   0.002;
3047       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
3048       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
3049     }
3050     else
3051     if (fParticleID==AliPID::kKaon)
3052     {
3053       //TOF kaons, 0.9 purity
3054       t = new TMatrixF(3,61);
3055       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
3056       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
3057       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
3058       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
3059       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.05;  (*t)(2,4)  =   0.05;
3060       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.05;  (*t)(2,5)  =   0.05;
3061       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.05;  (*t)(2,6)  =   0.05;
3062       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =   0.05;
3063       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.05;  (*t)(2,8)  =   0.05;
3064       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.05;  (*t)(2,9)  =   0.05;
3065       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.05;  (*t)(2,10)  =   0.05;
3066       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.026;  (*t)(2,11)  =   0.026;
3067       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.026;  (*t)(2,12)  =   0.026;
3068       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.026;  (*t)(2,13)  =   0.026;
3069       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.026;  (*t)(2,14)  =   0.026;
3070       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.026;  (*t)(2,15)  =   0.026;
3071       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.026;  (*t)(2,16)  =   0.026;
3072       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.024;  (*t)(2,17)  =   0.024;
3073       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.024;  (*t)(2,18)  =   0.024;
3074       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.024;  (*t)(2,19)  =   0.024;
3075       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.024;  (*t)(2,20)  =   0.024;
3076       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.024;  (*t)(2,21)  =   0.024;
3077       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.024;  (*t)(2,22)  =   0.022;
3078       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.020;
3079       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.026;  (*t)(2,24)  =   0.016;
3080       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.028;  (*t)(2,25)  =   0.014;
3081       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.028;  (*t)(2,26)  =   0.012;
3082       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.028;  (*t)(2,27)  =   0.010;
3083       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.028;  (*t)(2,28)  =   0.010;
3084       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.028;  (*t)(2,29)  =   0.008;
3085       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.028;  (*t)(2,30)  =   0.006;
3086       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.006;
3087       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.024;  (*t)(2,32)  =   0.004;
3088       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.020;  (*t)(2,33)  =   0.002;
3089       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.020;  (*t)(2,34)  =   0.002;
3090       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.018;  (*t)(2,35)  =   0.000;
3091       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.016;  (*t)(2,36)  =   0.000;
3092       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.014;  (*t)(2,37)  =   -0.002;
3093       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.014;  (*t)(2,38)  =   -0.004;
3094       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.012;  (*t)(2,39)  =   -0.004;
3095       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.010;  (*t)(2,40)  =   -0.006;
3096       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.000;  (*t)(2,41)  =   0.000;
3097       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.000;  (*t)(2,42)  =   0.000;
3098       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.000;  (*t)(2,43)  =   0.000;
3099       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.000;  (*t)(2,44)  =   0.000;
3100       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.000;  (*t)(2,45)  =   0.000;
3101       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.000;  (*t)(2,46)  =   0.000;
3102       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.000;  (*t)(2,47)  =   0.000;
3103       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.000;  (*t)(2,48)  =   0.000;
3104       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.000;  (*t)(2,49)  =   0.000;
3105       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.000;  (*t)(2,50)  =   0.000;
3106       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.000;  (*t)(2,51)  =   0.000;
3107       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.000;  (*t)(2,52)  =   0.000;
3108       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.000;  (*t)(2,53)  =   0.000;
3109       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
3110       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
3111       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
3112       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
3113       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
3114       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
3115       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
3116     }
3117     delete fTOFpidCuts;
3118     fTOFpidCuts=t;
3119   }
3120 }
3121
3122 //-----------------------------------------------------------------------
3123 //-----------------------------------------------------------------------
3124 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliAODTrack* track)
3125 {
3126   fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3127   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3128
3129   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3130
3131   if(! kTPC) return kFALSE;
3132
3133   fProbBayes = 0.0;
3134
3135 switch (fParticleID)
3136   {
3137     case AliPID::kPion:
3138       fProbBayes = probabilities[2];
3139       break;
3140     case AliPID::kKaon:
3141        fProbBayes = probabilities[3];
3142      break;
3143     case AliPID::kProton:
3144        fProbBayes = probabilities[4];
3145       break;
3146     case AliPID::kElectron:
3147        fProbBayes = probabilities[0];
3148      break;
3149     case AliPID::kMuon:
3150        fProbBayes = probabilities[1];
3151      break;
3152     case AliPID::kDeuteron:
3153        fProbBayes = probabilities[5];
3154      break;
3155     case AliPID::kTriton:
3156        fProbBayes = probabilities[6];
3157      break;
3158     case AliPID::kHe3:
3159        fProbBayes = probabilities[7];
3160      break;
3161    default:
3162       return kFALSE;
3163   }
3164
3165  if(fProbBayes > fParticleProbability){
3166     if(!fCutCharge)
3167       return kTRUE;
3168     else if (fCutCharge && fCharge * track->Charge() > 0)
3169       return kTRUE;
3170   }
3171   return kFALSE;
3172 }
3173 //-----------------------------------------------------------------------
3174 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
3175 {
3176   //cut on TPC bayesian pid
3177   //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3178
3179   //Bool_t statusMatchingHard = TPCTOFagree(track);
3180   //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3181   //     return kFALSE;
3182   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3183   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3184   //Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
3185
3186   if(! kTPC) return kFALSE;
3187
3188   //  Bool_t statusMatchingHard = 1;
3189   //  Float_t mismProb = 0;
3190   //  if(kTOF){
3191   //    statusMatchingHard = TPCTOFagree(track);
3192   //    mismProb = fBayesianResponse->GetTOFMismProb();
3193   //  }
3194   //  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3195   //       return kFALSE;
3196
3197   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3198
3199   fProbBayes = 0.0;
3200
3201   switch (fParticleID)
3202   {
3203     case AliPID::kPion:
3204       fProbBayes = probabilities[2];
3205       break;
3206     case AliPID::kKaon:
3207       fProbBayes = probabilities[3];
3208      break;
3209     case AliPID::kProton:
3210       fProbBayes = probabilities[4];
3211       break;
3212     case AliPID::kElectron:
3213        fProbBayes = probabilities[0];
3214      break;
3215     case AliPID::kMuon:
3216        fProbBayes = probabilities[1];
3217      break;
3218     case AliPID::kDeuteron:
3219        fProbBayes = probabilities[5];
3220      break;
3221     case AliPID::kTriton:
3222        fProbBayes = probabilities[6];
3223      break;
3224     case AliPID::kHe3:
3225        fProbBayes = probabilities[7];
3226      break;
3227     default:
3228       return kFALSE;
3229   }
3230
3231   if(fProbBayes > fParticleProbability)
3232     {
3233       if(!fCutCharge)
3234         return kTRUE;
3235       else if (fCutCharge && fCharge * track->GetSign() > 0)
3236         return kTRUE;
3237     }
3238   return kFALSE;
3239 }
3240 //-----------------------------------------------------------------------
3241 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliAODTrack* track)
3242 {
3243 //check is track passes bayesian combined TOF+TPC pid cut
3244   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3245                      (track->GetStatus() & AliESDtrack::kTIME) &&
3246                      (track->GetTOFsignal() > 12000) &&
3247                      (track->GetTOFsignal() < 100000);
3248
3249   if (! goodtrack)
3250        return kFALSE;
3251
3252   Bool_t statusMatchingHard = TPCTOFagree(track);
3253 //ciao
3254   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3255        return kFALSE;
3256
3257  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3258   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3259
3260   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3261
3262   fProbBayes = 0.0;
3263
3264  switch (fParticleID)
3265   {
3266     case AliPID::kPion:
3267       fProbBayes = probabilities[2];
3268       break;
3269     case AliPID::kKaon:
3270        fProbBayes = probabilities[3];
3271      break;
3272     case AliPID::kProton:
3273        fProbBayes = probabilities[4];
3274       break;
3275     case AliPID::kElectron:
3276        fProbBayes = probabilities[0];
3277      break;
3278     case AliPID::kMuon:
3279        fProbBayes = probabilities[1];
3280      break;
3281     case AliPID::kDeuteron:
3282        fProbBayes = probabilities[5];
3283      break;
3284     case AliPID::kTriton:
3285        fProbBayes = probabilities[6];
3286      break;
3287     case AliPID::kHe3:
3288        fProbBayes = probabilities[7];
3289      break;
3290    default:
3291       return kFALSE;
3292   }
3293
3294  if(fProbBayes > fParticleProbability && mismProb < 0.5){
3295     if(!fCutCharge)
3296       return kTRUE;
3297     else if (fCutCharge && fCharge * track->Charge() > 0)
3298       return kTRUE;
3299   }
3300   return kFALSE;
3301
3302 }
3303 //-----------------------------------------------------------------------
3304 // part added by F. Noferini (some methods)
3305 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
3306 {
3307   //check is track passes bayesian combined TOF+TPC pid cut
3308   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3309                      (track->GetStatus() & AliESDtrack::kTIME) &&
3310                      (track->GetTOFsignal() > 12000) &&
3311                      (track->GetTOFsignal() < 100000) &&
3312                      (track->GetIntegratedLength() > 365);
3313
3314   if (! goodtrack)
3315        return kFALSE;
3316
3317   if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3318
3319   Bool_t statusMatchingHard = TPCTOFagree(track);
3320   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3321        return kFALSE;
3322
3323   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3324   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3325
3326   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3327
3328   fProbBayes = 0.0;
3329
3330   switch (fParticleID)
3331   {
3332     case AliPID::kPion:
3333       fProbBayes = probabilities[2];
3334       break;
3335     case AliPID::kKaon:
3336        fProbBayes = probabilities[3];
3337      break;
3338     case AliPID::kProton:
3339        fProbBayes = probabilities[4];
3340       break;
3341     case AliPID::kElectron:
3342        fProbBayes = probabilities[0];
3343      break;
3344     case AliPID::kMuon:
3345        fProbBayes = probabilities[1];
3346      break;
3347     case AliPID::kDeuteron:
3348        fProbBayes = probabilities[5];
3349      break;
3350     case AliPID::kTriton:
3351        fProbBayes = probabilities[6];
3352      break;
3353     case AliPID::kHe3:
3354        fProbBayes = probabilities[7];
3355      break;
3356    default:
3357       return kFALSE;
3358   }
3359
3360   //  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);
3361   if(fProbBayes > fParticleProbability && mismProb < 0.5){
3362     if(!fCutCharge)
3363       return kTRUE;
3364     else if (fCutCharge && fCharge * track->GetSign() > 0)
3365       return kTRUE;
3366   }
3367   return kFALSE;
3368 }
3369
3370
3371 //-----------------------------------------------------------------------
3372  // part added by Natasha
3373 Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
3374 {
3375   //pid selection for heavy nuclei
3376   Bool_t select=kFALSE;
3377
3378   //if (!track) continue;
3379
3380   if (!track->GetInnerParam())
3381     return kFALSE;    //break;
3382
3383   const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
3384
3385   Double_t ptotTPC = tpcTrack->GetP();
3386   Double_t sigTPC = track->GetTPCsignal();
3387   Double_t dEdxBBA = 0.;
3388   Double_t dSigma = 0.;
3389
3390   switch (fParticleID)
3391   {
3392     case AliPID::kDeuteron:
3393       //pid=10;
3394       dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
3395           4.60e+00,
3396           8.9684e+00,
3397           1.640e-05,
3398           2.35e+00,
3399           2.35e+00);
3400       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
3401
3402       if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) )  )
3403       {select=kTRUE;}
3404       break;
3405
3406     case AliPID::kTriton:
3407       //pid=11;
3408       select=kFALSE;
3409       break;
3410
3411     case AliPID::kHe3:
3412       //pid=12;
3413       // ----- Pass 2 -------
3414       dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
3415           1.74962,
3416           27.4992,
3417           4.00313e-15,
3418           2.42485,
3419           8.31768);
3420       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
3421       if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
3422       {select=kTRUE;}
3423       break;
3424
3425     case AliPID::kAlpha:
3426       //pid=13;
3427       select=kFALSE;
3428       break;
3429
3430     default:
3431       return kFALSE;
3432   }
3433
3434   return select;
3435 }
3436 // end part added by Natasha
3437 //-----------------------------------------------------------------------
3438 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
3439  //set priors for the bayesian pid selection
3440   fCurrCentr = centrCur;
3441
3442   fBinLimitPID[0] = 0.300000;
3443   fBinLimitPID[1] = 0.400000;
3444   fBinLimitPID[2] = 0.500000;
3445   fBinLimitPID[3] = 0.600000;
3446   fBinLimitPID[4] = 0.700000;
3447   fBinLimitPID[5] = 0.800000;
3448   fBinLimitPID[6] = 0.900000;
3449   fBinLimitPID[7] = 1.000000;
3450   fBinLimitPID[8] = 1.200000;
3451   fBinLimitPID[9] = 1.400000;
3452   fBinLimitPID[10] = 1.600000;
3453   fBinLimitPID[11] = 1.800000;
3454   fBinLimitPID[12] = 2.000000;
3455   fBinLimitPID[13] = 2.200000;
3456   fBinLimitPID[14] = 2.400000;
3457   fBinLimitPID[15] = 2.600000;
3458   fBinLimitPID[16] = 2.800000;
3459   fBinLimitPID[17] = 3.000000;
3460
3461   // 0-10%
3462   if(centrCur < 10){
3463       fC[0][0] = 0.005;
3464       fC[0][1] = 0.005;
3465       fC[0][2] = 1.0000;
3466       fC[0][3] = 0.0010;
3467       fC[0][4] = 0.0010;
3468
3469       fC[1][0] = 0.005;
3470       fC[1][1] = 0.005;
3471       fC[1][2] = 1.0000;
3472       fC[1][3] = 0.0168;
3473       fC[1][4] = 0.0122;
3474
3475       fC[2][0] = 0.005;
3476       fC[2][1] = 0.005;
3477       fC[2][2] = 1.0000;
3478       fC[2][3] = 0.0272;
3479       fC[2][4] = 0.0070;
3480
3481       fC[3][0] = 0.005;
3482       fC[3][1] = 0.005;
3483       fC[3][2] = 1.0000;
3484       fC[3][3] = 0.0562;
3485       fC[3][4] = 0.0258;
3486
3487       fC[4][0] = 0.005;
3488       fC[4][1] = 0.005;
3489       fC[4][2] = 1.0000;
3490       fC[4][3] = 0.0861;
3491       fC[4][4] = 0.0496;
3492
3493       fC[5][0] = 0.005;
3494       fC[5][1] = 0.005;
3495       fC[5][2] = 1.0000;
3496       fC[5][3] = 0.1168;
3497       fC[5][4] = 0.0740;
3498
3499       fC[6][0] = 0.005;
3500       fC[6][1] = 0.005;
3501       fC[6][2] = 1.0000;
3502       fC[6][3] = 0.1476;
3503       fC[6][4] = 0.0998;
3504
3505       fC[7][0] = 0.005;
3506       fC[7][1] = 0.005;
3507       fC[7][2] = 1.0000;
3508       fC[7][3] = 0.1810;
3509       fC[7][4] = 0.1296;
3510
3511       fC[8][0] = 0.005;
3512       fC[8][1] = 0.005;
3513       fC[8][2] = 1.0000;
3514       fC[8][3] = 0.2240;
3515       fC[8][4] = 0.1827;
3516
3517       fC[9][0] = 0.005;
3518       fC[9][1] = 0.005;
3519       fC[9][2] = 1.0000;
3520       fC[9][3] = 0.2812;
3521       fC[9][4] = 0.2699;
3522
3523       fC[10][0] = 0.005;
3524       fC[10][1] = 0.005;
3525       fC[10][2] = 1.0000;
3526       fC[10][3] = 0.3328;
3527       fC[10][4] = 0.3714;
3528
3529       fC[11][0] = 0.005;
3530       fC[11][1] = 0.005;
3531       fC[11][2] = 1.0000;
3532       fC[11][3] = 0.3780;
3533       fC[11][4] = 0.4810;
3534
3535       fC[12][0] = 0.005;
3536       fC[12][1] = 0.005;
3537       fC[12][2] = 1.0000;
3538       fC[12][3] = 0.4125;
3539       fC[12][4] = 0.5771;
3540
3541       fC[13][0] = 0.005;
3542       fC[13][1] = 0.005;
3543       fC[13][2] = 1.0000;
3544       fC[13][3] = 0.4486;
3545       fC[13][4] = 0.6799;
3546
3547       fC[14][0] = 0.005;
3548       fC[14][1] = 0.005;
3549       fC[14][2] = 1.0000;
3550       fC[14][3] = 0.4840;
3551       fC[14][4] = 0.7668;
3552
3553       fC[15][0] = 0.005;
3554       fC[15][1] = 0.005;
3555       fC[15][2] = 1.0000;
3556       fC[15][3] = 0.4971;
3557       fC[15][4] = 0.8288;
3558
3559       fC[16][0] = 0.005;
3560       fC[16][1] = 0.005;
3561       fC[16][2] = 1.0000;
3562       fC[16][3] = 0.4956;
3563       fC[16][4] = 0.8653;
3564
3565       fC[17][0] = 0.005;
3566       fC[17][1] = 0.005;
3567       fC[17][2] = 1.0000;
3568       fC[17][3] = 0.5173;
3569       fC[17][4] = 0.9059;
3570   }
3571   // 10-20%
3572   else if(centrCur < 20){
3573      fC[0][0] = 0.005;
3574       fC[0][1] = 0.005;
3575       fC[0][2] = 1.0000;
3576       fC[0][3] = 0.0010;
3577       fC[0][4] = 0.0010;
3578
3579       fC[1][0] = 0.005;
3580       fC[1][1] = 0.005;
3581       fC[1][2] = 1.0000;
3582       fC[1][3] = 0.0132;
3583       fC[1][4] = 0.0088;
3584
3585       fC[2][0] = 0.005;
3586       fC[2][1] = 0.005;
3587       fC[2][2] = 1.0000;
3588       fC[2][3] = 0.0283;
3589       fC[2][4] = 0.0068;
3590
3591       fC[3][0] = 0.005;
3592       fC[3][1] = 0.005;
3593       fC[3][2] = 1.0000;
3594       fC[3][3] = 0.0577;
3595       fC[3][4] = 0.0279;
3596
3597       fC[4][0] = 0.005;
3598       fC[4][1] = 0.005;
3599       fC[4][2] = 1.0000;
3600       fC[4][3] = 0.0884;
3601       fC[4][4] = 0.0534;
3602
3603       fC[5][0] = 0.005;
3604       fC[5][1] = 0.005;
3605       fC[5][2] = 1.0000;
3606       fC[5][3] = 0.1179;
3607       fC[5][4] = 0.0794;
3608
3609       fC[6][0] = 0.005;
3610       fC[6][1] = 0.005;
3611       fC[6][2] = 1.0000;
3612       fC[6][3] = 0.1480;
3613       fC[6][4] = 0.1058;
3614
3615       fC[7][0] = 0.005;
3616       fC[7][1] = 0.005;
3617       fC[7][2] = 1.0000;
3618       fC[7][3] = 0.1807;
3619       fC[7][4] = 0.1366;
3620
3621       fC[8][0] = 0.005;
3622       fC[8][1] = 0.005;
3623       fC[8][2] = 1.0000;
3624       fC[8][3] = 0.2219;
3625       fC[8][4] = 0.1891;
3626
3627       fC[9][0] = 0.005;
3628       fC[9][1] = 0.005;
3629       fC[9][2] = 1.0000;
3630       fC[9][3] = 0.2804;
3631       fC[9][4] = 0.2730;
3632
3633       fC[10][0] = 0.005;
3634       fC[10][1] = 0.005;
3635       fC[10][2] = 1.0000;
3636       fC[10][3] = 0.3283;
3637       fC[10][4] = 0.3660;
3638
3639       fC[11][0] = 0.005;
3640       fC[11][1] = 0.005;
3641       fC[11][2] = 1.0000;
3642       fC[11][3] = 0.3710;
3643       fC[11][4] = 0.4647;
3644
3645       fC[12][0] = 0.005;
3646       fC[12][1] = 0.005;
3647       fC[12][2] = 1.0000;
3648       fC[12][3] = 0.4093;
3649       fC[12][4] = 0.5566;
3650
3651       fC[13][0] = 0.005;
3652       fC[13][1] = 0.005;
3653       fC[13][2] = 1.0000;
3654       fC[13][3] = 0.4302;
3655       fC[13][4] = 0.6410;
3656
3657       fC[14][0] = 0.005;
3658       fC[14][1] = 0.005;
3659       fC[14][2] = 1.0000;
3660       fC[14][3] = 0.4649;
3661       fC[14][4] = 0.7055;
3662
3663       fC[15][0] = 0.005;
3664       fC[15][1] = 0.005;
3665       fC[15][2] = 1.0000;
3666       fC[15][3] = 0.4523;
3667       fC[15][4] = 0.7440;
3668
3669       fC[16][0] = 0.005;
3670       fC[16][1] = 0.005;
3671       fC[16][2] = 1.0000;
3672       fC[16][3] = 0.4591;
3673       fC[16][4] = 0.7799;
3674
3675       fC[17][0] = 0.005;
3676       fC[17][1] = 0.005;
3677       fC[17][2] = 1.0000;
3678       fC[17][3] = 0.4804;
3679       fC[17][4] = 0.8218;
3680   }
3681   // 20-30%
3682   else if(centrCur < 30){
3683      fC[0][0] = 0.005;
3684       fC[0][1] = 0.005;
3685       fC[0][2] = 1.0000;
3686       fC[0][3] = 0.0010;
3687       fC[0][4] = 0.0010;
3688
3689       fC[1][0] = 0.005;
3690       fC[1][1] = 0.005;
3691       fC[1][2] = 1.0000;
3692       fC[1][3] = 0.0102;
3693       fC[1][4] = 0.0064;
3694
3695       fC[2][0] = 0.005;
3696       fC[2][1] = 0.005;
3697       fC[2][2] = 1.0000;
3698       fC[2][3] = 0.0292;
3699       fC[2][4] = 0.0066;
3700
3701       fC[3][0] = 0.005;
3702       fC[3][1] = 0.005;
3703       fC[3][2] = 1.0000;
3704       fC[3][3] = 0.0597;
3705       fC[3][4] = 0.0296;
3706
3707       fC[4][0] = 0.005;
3708       fC[4][1] = 0.005;
3709       fC[4][2] = 1.0000;
3710       fC[4][3] = 0.0900;
3711       fC[4][4] = 0.0589;
3712
3713       fC[5][0] = 0.005;
3714       fC[5][1] = 0.005;
3715       fC[5][2] = 1.0000;
3716       fC[5][3] = 0.1199;
3717       fC[5][4] = 0.0859;
3718
3719       fC[6][0] = 0.005;
3720       fC[6][1] = 0.005;
3721       fC[6][2] = 1.0000;
3722       fC[6][3] = 0.1505;
3723       fC[6][4] = 0.1141;
3724
3725       fC[7][0] = 0.005;
3726       fC[7][1] = 0.005;
3727       fC[7][2] = 1.0000;
3728       fC[7][3] = 0.1805;
3729       fC[7][4] = 0.1454;
3730
3731       fC[8][0] = 0.005;
3732       fC[8][1] = 0.005;
3733       fC[8][2] = 1.0000;
3734       fC[8][3] = 0.2221;
3735       fC[8][4] = 0.2004;
3736
3737       fC[9][0] = 0.005;
3738       fC[9][1] = 0.005;
3739       fC[9][2] = 1.0000;
3740       fC[9][3] = 0.2796;
3741       fC[9][4] = 0.2838;
3742
3743       fC[10][0] = 0.005;
3744       fC[10][1] = 0.005;
3745       fC[10][2] = 1.0000;
3746       fC[10][3] = 0.3271;
3747       fC[10][4] = 0.3682;
3748
3749       fC[11][0] = 0.005;
3750       fC[11][1] = 0.005;
3751       fC[11][2] = 1.0000;
3752       fC[11][3] = 0.3648;
3753       fC[11][4] = 0.4509;
3754
3755       fC[12][0] = 0.005;
3756       fC[12][1] = 0.005;
3757       fC[12][2] = 1.0000;
3758       fC[12][3] = 0.3988;
3759       fC[12][4] = 0.5339;
3760
3761       fC[13][0] = 0.005;
3762       fC[13][1] = 0.005;
3763       fC[13][2] = 1.0000;
3764       fC[13][3] = 0.4315;
3765       fC[13][4] = 0.5995;
3766
3767       fC[14][0] = 0.005;
3768       fC[14][1] = 0.005;
3769       fC[14][2] = 1.0000;
3770       fC[14][3] = 0.4548;
3771       fC[14][4] = 0.6612;
3772
3773       fC[15][0] = 0.005;
3774       fC[15][1] = 0.005;
3775       fC[15][2] = 1.0000;
3776       fC[15][3] = 0.4744;
3777       fC[15][4] = 0.7060;
3778
3779       fC[16][0] = 0.005;
3780       fC[16][1] = 0.005;
3781       fC[16][2] = 1.0000;
3782       fC[16][3] = 0.4899;
3783       fC[16][4] = 0.7388;
3784
3785       fC[17][0] = 0.005;
3786       fC[17][1] = 0.005;
3787       fC[17][2] = 1.0000;
3788       fC[17][3] = 0.4411;
3789       fC[17][4] = 0.7293;
3790   }
3791   // 30-40%
3792   else if(centrCur < 40){
3793       fC[0][0] = 0.005;
3794       fC[0][1] = 0.005;
3795       fC[0][2] = 1.0000;
3796       fC[0][3] = 0.0010;
3797       fC[0][4] = 0.0010;
3798
3799       fC[1][0] = 0.005;
3800       fC[1][1] = 0.005;
3801       fC[1][2] = 1.0000;
3802       fC[1][3] = 0.0102;
3803       fC[1][4] = 0.0048;
3804
3805       fC[2][0] = 0.005;
3806       fC[2][1] = 0.005;
3807       fC[2][2] = 1.0000;
3808       fC[2][3] = 0.0306;
3809       fC[2][4] = 0.0079;
3810
3811       fC[3][0] = 0.005;
3812       fC[3][1] = 0.005;
3813       fC[3][2] = 1.0000;
3814       fC[3][3] = 0.0617;
3815       fC[3][4] = 0.0338;
3816
3817       fC[4][0] = 0.005;
3818       fC[4][1] = 0.005;
3819       fC[4][2] = 1.0000;
3820       fC[4][3] = 0.0920;
3821       fC[4][4] = 0.0652;
3822
3823       fC[5][0] = 0.005;
3824       fC[5][1] = 0.005;
3825       fC[5][2] = 1.0000;
3826       fC[5][3] = 0.1211;
3827       fC[5][4] = 0.0955;
3828
3829       fC[6][0] = 0.005;
3830       fC[6][1] = 0.005;
3831       fC[6][2] = 1.0000;
3832       fC[6][3] = 0.1496;
3833       fC[6][4] = 0.1242;
3834
3835       fC[7][0] = 0.005;
3836       fC[7][1] = 0.005;
3837       fC[7][2] = 1.0000;
3838       fC[7][3] = 0.1807;
3839       fC[7][4] = 0.1576;
3840
3841       fC[8][0] = 0.005;
3842       fC[8][1] = 0.005;
3843       fC[8][2] = 1.0000;
3844       fC[8][3] = 0.2195;
3845       fC[8][4] = 0.2097;
3846
3847       fC[9][0] = 0.005;
3848       fC[9][1] = 0.005;
3849       fC[9][2] = 1.0000;
3850       fC[9][3] = 0.2732;
3851       fC[9][4] = 0.2884;
3852
3853       fC[10][0] = 0.005;
3854       fC[10][1] = 0.005;
3855       fC[10][2] = 1.0000;
3856       fC[10][3] = 0.3204;
3857       fC[10][4] = 0.3679;
3858
3859       fC[11][0] = 0.005;
3860       fC[11][1] = 0.005;
3861       fC[11][2] = 1.0000;
3862       fC[11][3] = 0.3564;
3863       fC[11][4] = 0.4449;
3864
3865       fC[12][0] = 0.005;
3866       fC[12][1] = 0.005;
3867       fC[12][2] = 1.0000;
3868       fC[12][3] = 0.3791;
3869       fC[12][4] = 0.5052;
3870
3871       fC[13][0] = 0.005;
3872       fC[13][1] = 0.005;
3873       fC[13][2] = 1.0000;
3874       fC[13][3] = 0.4062;
3875       fC[13][4] = 0.5647;
3876
3877       fC[14][0] = 0.005;
3878       fC[14][1] = 0.005;
3879       fC[14][2] = 1.0000;
3880       fC[14][3] = 0.4234;
3881       fC[14][4] = 0.6203;
3882
3883       fC[15][0] = 0.005;
3884       fC[15][1] = 0.005;
3885       fC[15][2] = 1.0000;
3886       fC[15][3] = 0.4441;
3887       fC[15][4] = 0.6381;
3888
3889       fC[16][0] = 0.005;
3890       fC[16][1] = 0.005;
3891       fC[16][2] = 1.0000;
3892       fC[16][3] = 0.4629;
3893       fC[16][4] = 0.6496;
3894
3895       fC[17][0] = 0.005;
3896       fC[17][1] = 0.005;
3897       fC[17][2] = 1.0000;
3898       fC[17][3] = 0.4293;
3899       fC[17][4] = 0.6491;
3900   }
3901   // 40-50%
3902   else if(centrCur < 50){
3903       fC[0][0] = 0.005;
3904       fC[0][1] = 0.005;
3905       fC[0][2] = 1.0000;
3906       fC[0][3] = 0.0010;
3907       fC[0][4] = 0.0010;
3908
3909       fC[1][0] = 0.005;
3910       fC[1][1] = 0.005;
3911       fC[1][2] = 1.0000;
3912       fC[1][3] = 0.0093;
3913       fC[1][4] = 0.0057;
3914
3915       fC[2][0] = 0.005;
3916       fC[2][1] = 0.005;
3917       fC[2][2] = 1.0000;
3918       fC[2][3] = 0.0319;
3919       fC[2][4] = 0.0075;
3920
3921       fC[3][0] = 0.005;
3922       fC[3][1] = 0.005;
3923       fC[3][2] = 1.0000;
3924       fC[3][3] = 0.0639;
3925       fC[3][4] = 0.0371;
3926
3927       fC[4][0] = 0.005;
3928       fC[4][1] = 0.005;
3929       fC[4][2] = 1.0000;
3930       fC[4][3] = 0.0939;
3931       fC[4][4] = 0.0725;
3932
3933       fC[5][0] = 0.005;
3934       fC[5][1] = 0.005;
3935       fC[5][2] = 1.0000;
3936       fC[5][3] = 0.1224;
3937       fC[5][4] = 0.1045;
3938
3939       fC[6][0] = 0.005;
3940       fC[6][1] = 0.005;
3941       fC[6][2] = 1.0000;
3942       fC[6][3] = 0.1520;
3943       fC[6][4] = 0.1387;
3944
3945       fC[7][0] = 0.005;
3946       fC[7][1] = 0.005;
3947       fC[7][2] = 1.0000;
3948       fC[7][3] = 0.1783;
3949       fC[7][4] = 0.1711;
3950
3951       fC[8][0] = 0.005;
3952       fC[8][1] = 0.005;
3953       fC[8][2] = 1.0000;
3954       fC[8][3] = 0.2202;
3955       fC[8][4] = 0.2269;
3956
3957       fC[9][0] = 0.005;
3958       fC[9][1] = 0.005;
3959       fC[9][2] = 1.0000;
3960       fC[9][3] = 0.2672;
3961       fC[9][4] = 0.2955;
3962
3963       fC[10][0] = 0.005;
3964       fC[10][1] = 0.005;
3965       fC[10][2] = 1.0000;
3966       fC[10][3] = 0.3191;
3967       fC[10][4] = 0.3676;
3968
3969       fC[11][0] = 0.005;
3970       fC[11][1] = 0.005;
3971       fC[11][2] = 1.0000;
3972       fC[11][3] = 0.3434;
3973       fC[11][4] = 0.4321;
3974
3975       fC[12][0] = 0.005;
3976       fC[12][1] = 0.005;
3977       fC[12][2] = 1.0000;
3978       fC[12][3] = 0.3692;
3979       fC[12][4] = 0.4879;
3980
3981       fC[13][0] = 0.005;
3982       fC[13][1] = 0.005;
3983       fC[13][2] = 1.0000;
3984       fC[13][3] = 0.3993;
3985       fC[13][4] = 0.5377;
3986
3987       fC[14][0] = 0.005;
3988       fC[14][1] = 0.005;
3989       fC[14][2] = 1.0000;
3990       fC[14][3] = 0.3818;
3991       fC[14][4] = 0.5547;
3992
3993       fC[15][0] = 0.005;
3994       fC[15][1] = 0.005;
3995       fC[15][2] = 1.0000;
3996       fC[15][3] = 0.4003;
3997       fC[15][4] = 0.5484;
3998
3999       fC[16][0] = 0.005;
4000       fC[16][1] = 0.005;
4001       fC[16][2] = 1.0000;
4002       fC[16][3] = 0.4281;
4003       fC[16][4] = 0.5383;
4004
4005       fC[17][0] = 0.005;
4006       fC[17][1] = 0.005;
4007       fC[17][2] = 1.0000;
4008       fC[17][3] = 0.3960;
4009       fC[17][4] = 0.5374;
4010   }
4011   // 50-60%
4012   else if(centrCur < 60){
4013       fC[0][0] = 0.005;
4014       fC[0][1] = 0.005;
4015       fC[0][2] = 1.0000;
4016       fC[0][3] = 0.0010;
4017       fC[0][4] = 0.0010;
4018
4019       fC[1][0] = 0.005;
4020       fC[1][1] = 0.005;
4021       fC[1][2] = 1.0000;
4022       fC[1][3] = 0.0076;
4023       fC[1][4] = 0.0032;
4024
4025       fC[2][0] = 0.005;
4026       fC[2][1] = 0.005;
4027       fC[2][2] = 1.0000;
4028       fC[2][3] = 0.0329;
4029       fC[2][4] = 0.0085;
4030
4031       fC[3][0] = 0.005;
4032       fC[3][1] = 0.005;
4033       fC[3][2] = 1.0000;
4034       fC[3][3] = 0.0653;
4035       fC[3][4] = 0.0423;
4036
4037       fC[4][0] = 0.005;
4038       fC[4][1] = 0.005;
4039       fC[4][2] = 1.0000;
4040       fC[4][3] = 0.0923;
4041       fC[4][4] = 0.0813;
4042
4043       fC[5][0] = 0.005;
4044       fC[5][1] = 0.005;
4045       fC[5][2] = 1.0000;
4046       fC[5][3] = 0.1219;
4047       fC[5][4] = 0.1161;
4048
4049       fC[6][0] = 0.005;
4050       fC[6][1] = 0.005;
4051       fC[6][2] = 1.0000;
4052       fC[6][3] = 0.1519;
4053       fC[6][4] = 0.1520;
4054
4055       fC[7][0] = 0.005;
4056       fC[7][1] = 0.005;
4057       fC[7][2] = 1.0000;
4058       fC[7][3] = 0.1763;
4059       fC[7][4] = 0.1858;
4060
4061       fC[8][0] = 0.005;
4062       fC[8][1] = 0.005;
4063       fC[8][2] = 1.0000;
4064       fC[8][3] = 0.2178;
4065       fC[8][4] = 0.2385;
4066
4067       fC[9][0] = 0.005;
4068       fC[9][1] = 0.005;
4069       fC[9][2] = 1.0000;
4070       fC[9][3] = 0.2618;
4071       fC[9][4] = 0.3070;
4072
4073       fC[10][0] = 0.005;
4074       fC[10][1] = 0.005;
4075       fC[10][2] = 1.0000;
4076       fC[10][3] = 0.3067;
4077       fC[10][4] = 0.3625;
4078
4079       fC[11][0] = 0.005;
4080       fC[11][1] = 0.005;
4081       fC[11][2] = 1.0000;
4082       fC[11][3] = 0.3336;
4083       fC[11][4] = 0.4188;
4084
4085       fC[12][0] = 0.005;
4086       fC[12][1] = 0.005;
4087       fC[12][2] = 1.0000;
4088       fC[12][3] = 0.3706;
4089       fC[12][4] = 0.4511;
4090
4091       fC[13][0] = 0.005;
4092       fC[13][1] = 0.005;
4093       fC[13][2] = 1.0000;
4094       fC[13][3] = 0.3765;
4095       fC[13][4] = 0.4729;
4096
4097       fC[14][0] = 0.005;
4098       fC[14][1] = 0.005;
4099       fC[14][2] = 1.0000;
4100       fC[14][3] = 0.3942;
4101       fC[14][4] = 0.4855;
4102
4103       fC[15][0] = 0.005;
4104       fC[15][1] = 0.005;
4105       fC[15][2] = 1.0000;
4106       fC[15][3] = 0.4051;
4107       fC[15][4] = 0.4762;
4108
4109       fC[16][0] = 0.005;
4110       fC[16][1] = 0.005;
4111       fC[16][2] = 1.0000;
4112       fC[16][3] = 0.3843;
4113       fC[16][4] = 0.4763;
4114
4115       fC[17][0] = 0.005;
4116       fC[17][1] = 0.005;
4117       fC[17][2] = 1.0000;
4118       fC[17][3] = 0.4237;
4119       fC[17][4] = 0.4773;
4120   }
4121   // 60-70%
4122   else if(centrCur < 70){
4123          fC[0][0] = 0.005;
4124       fC[0][1] = 0.005;
4125       fC[0][2] = 1.0000;
4126       fC[0][3] = 0.0010;
4127       fC[0][4] = 0.0010;
4128
4129       fC[1][0] = 0.005;
4130       fC[1][1] = 0.005;
4131       fC[1][2] = 1.0000;
4132       fC[1][3] = 0.0071;
4133       fC[1][4] = 0.0012;
4134
4135       fC[2][0] = 0.005;
4136       fC[2][1] = 0.005;
4137       fC[2][2] = 1.0000;
4138       fC[2][3] = 0.0336;
4139       fC[2][4] = 0.0097;
4140
4141       fC[3][0] = 0.005;
4142       fC[3][1] = 0.005;
4143       fC[3][2] = 1.0000;
4144       fC[3][3] = 0.0662;
4145       fC[3][4] = 0.0460;
4146
4147       fC[4][0] = 0.005;
4148       fC[4][1] = 0.005;
4149       fC[4][2] = 1.0000;
4150       fC[4][3] = 0.0954;
4151       fC[4][4] = 0.0902;
4152
4153       fC[5][0] = 0.005;
4154       fC[5][1] = 0.005;
4155       fC[5][2] = 1.0000;
4156       fC[5][3] = 0.1181;
4157       fC[5][4] = 0.1306;
4158
4159       fC[6][0] = 0.005;
4160       fC[6][1] = 0.005;
4161       fC[6][2] = 1.0000;
4162       fC[6][3] = 0.1481;
4163       fC[6][4] = 0.1662;
4164
4165       fC[7][0] = 0.005;
4166       fC[7][1] = 0.005;
4167       fC[7][2] = 1.0000;
4168       fC[7][3] = 0.1765;
4169       fC[7][4] = 0.1963;
4170
4171       fC[8][0] = 0.005;
4172       fC[8][1] = 0.005;
4173       fC[8][2] = 1.0000;
4174       fC[8][3] = 0.2155;
4175       fC[8][4] = 0.2433;
4176
4177       fC[9][0] = 0.005;
4178       fC[9][1] = 0.005;
4179       fC[9][2] = 1.0000;
4180       fC[9][3] = 0.2580;
4181       fC[9][4] = 0.3022;
4182
4183       fC[10][0] = 0.005;
4184       fC[10][1] = 0.005;
4185       fC[10][2] = 1.0000;
4186       fC[10][3] = 0.2872;
4187       fC[10][4] = 0.3481;
4188
4189       fC[11][0] = 0.005;
4190       fC[11][1] = 0.005;
4191       fC[11][2] = 1.0000;
4192       fC[11][3] = 0.3170;
4193       fC[11][4] = 0.3847;
4194
4195       fC[12][0] = 0.005;
4196       fC[12][1] = 0.005;
4197       fC[12][2] = 1.0000;
4198       fC[12][3] = 0.3454;
4199       fC[12][4] = 0.4258;
4200
4201       fC[13][0] = 0.005;
4202       fC[13][1] = 0.005;
4203       fC[13][2] = 1.0000;
4204       fC[13][3] = 0.3580;
4205       fC[13][4] = 0.4299;
4206
4207       fC[14][0] = 0.005;
4208       fC[14][1] = 0.005;
4209       fC[14][2] = 1.0000;
4210       fC[14][3] = 0.3903;
4211       fC[14][4] = 0.4326;
4212
4213       fC[15][0] = 0.005;
4214       fC[15][1] = 0.005;
4215       fC[15][2] = 1.0000;
4216       fC[15][3] = 0.3690;
4217       fC[15][4] = 0.4491;
4218
4219       fC[16][0] = 0.005;
4220       fC[16][1] = 0.005;
4221       fC[16][2] = 1.0000;
4222       fC[16][3] = 0.4716;
4223       fC[16][4] = 0.4298;
4224
4225       fC[17][0] = 0.005;
4226       fC[17][1] = 0.005;
4227       fC[17][2] = 1.0000;
4228       fC[17][3] = 0.3875;
4229       fC[17][4] = 0.4083;
4230   }
4231   // 70-80%
4232   else if(centrCur < 80){
4233       fC[0][0] = 0.005;
4234       fC[0][1] = 0.005;
4235       fC[0][2] = 1.0000;
4236       fC[0][3] = 0.0010;
4237       fC[0][4] = 0.0010;
4238
4239       fC[1][0] = 0.005;
4240       fC[1][1] = 0.005;
4241       fC[1][2] = 1.0000;
4242       fC[1][3] = 0.0075;
4243       fC[1][4] = 0.0007;
4244
4245       fC[2][0] = 0.005;
4246       fC[2][1] = 0.005;
4247       fC[2][2] = 1.0000;
4248       fC[2][3] = 0.0313;
4249       fC[2][4] = 0.0124;
4250
4251       fC[3][0] = 0.005;
4252       fC[3][1] = 0.005;
4253       fC[3][2] = 1.0000;
4254       fC[3][3] = 0.0640;
4255       fC[3][4] = 0.0539;
4256
4257       fC[4][0] = 0.005;
4258       fC[4][1] = 0.005;
4259       fC[4][2] = 1.0000;
4260       fC[4][3] = 0.0923;
4261       fC[4][4] = 0.0992;
4262
4263       fC[5][0] = 0.005;
4264       fC[5][1] = 0.005;
4265       fC[5][2] = 1.0000;
4266       fC[5][3] = 0.1202;
4267       fC[5][4] = 0.1417;
4268
4269       fC[6][0] = 0.005;
4270       fC[6][1] = 0.005;
4271       fC[6][2] = 1.0000;
4272       fC[6][3] = 0.1413;
4273       fC[6][4] = 0.1729;
4274
4275       fC[7][0] = 0.005;
4276       fC[7][1] = 0.005;
4277       fC[7][2] = 1.0000;
4278       fC[7][3] = 0.1705;
4279       fC[7][4] = 0.1999;
4280
4281       fC[8][0] = 0.005;
4282       fC[8][1] = 0.005;
4283       fC[8][2] = 1.0000;
4284       fC[8][3] = 0.2103;
4285       fC[8][4] = 0.2472;
4286
4287       fC[9][0] = 0.005;
4288       fC[9][1] = 0.005;
4289       fC[9][2] = 1.0000;
4290       fC[9][3] = 0.2373;
4291       fC[9][4] = 0.2916;
4292
4293       fC[10][0] = 0.005;
4294       fC[10][1] = 0.005;
4295       fC[10][2] = 1.0000;
4296       fC[10][3] = 0.2824;
4297       fC[10][4] = 0.3323;
4298
4299       fC[11][0] = 0.005;
4300       fC[11][1] = 0.005;
4301       fC[11][2] = 1.0000;
4302       fC[11][3] = 0.3046;
4303       fC[11][4] = 0.3576;
4304
4305       fC[12][0] = 0.005;
4306       fC[12][1] = 0.005;
4307       fC[12][2] = 1.0000;
4308       fC[12][3] = 0.3585;
4309       fC[12][4] = 0.4003;
4310
4311       fC[13][0] = 0.005;
4312       fC[13][1] = 0.005;
4313       fC[13][2] = 1.0000;
4314       fC[13][3] = 0.3461;
4315       fC[13][4] = 0.3982;
4316
4317       fC[14][0] = 0.005;
4318       fC[14][1] = 0.005;
4319       fC[14][2] = 1.0000;
4320       fC[14][3] = 0.3362;
4321       fC[14][4] = 0.3776;
4322
4323       fC[15][0] = 0.005;
4324       fC[15][1] = 0.005;
4325       fC[15][2] = 1.0000;
4326       fC[15][3] = 0.3071;
4327       fC[15][4] = 0.3500;
4328
4329       fC[16][0] = 0.005;
4330       fC[16][1] = 0.005;
4331       fC[16][2] = 1.0000;
4332       fC[16][3] = 0.2914;
4333       fC[16][4] = 0.3937;
4334
4335       fC[17][0] = 0.005;
4336       fC[17][1] = 0.005;
4337       fC[17][2] = 1.0000;
4338       fC[17][3] = 0.3727;
4339       fC[17][4] = 0.3877;
4340   }
4341   // 80-100%
4342   else{
4343       fC[0][0] = 0.005;
4344       fC[0][1] = 0.005;
4345       fC[0][2] = 1.0000;
4346       fC[0][3] = 0.0010;
4347       fC[0][4] = 0.0010;
4348
4349       fC[1][0] = 0.005;
4350       fC[1][1] = 0.005;
4351       fC[1][2] = 1.0000;
4352       fC[1][3] = 0.0060;
4353       fC[1][4] = 0.0035;
4354
4355       fC[2][0] = 0.005;
4356       fC[2][1] = 0.005;
4357       fC[2][2] = 1.0000;
4358       fC[2][3] = 0.0323;
4359       fC[2][4] = 0.0113;
4360
4361       fC[3][0] = 0.005;
4362       fC[3][1] = 0.005;
4363       fC[3][2] = 1.0000;
4364       fC[3][3] = 0.0609;
4365       fC[3][4] = 0.0653;
4366
4367       fC[4][0] = 0.005;
4368       fC[4][1] = 0.005;
4369       fC[4][2] = 1.0000;
4370       fC[4][3] = 0.0922;
4371       fC[4][4] = 0.1076;
4372
4373       fC[5][0] = 0.005;
4374       fC[5][1] = 0.005;
4375       fC[5][2] = 1.0000;
4376       fC[5][3] = 0.1096;
4377       fC[5][4] = 0.1328;
4378
4379       fC[6][0] = 0.005;
4380       fC[6][1] = 0.005;
4381       fC[6][2] = 1.0000;
4382       fC[6][3] = 0.1495;
4383       fC[6][4] = 0.1779;
4384
4385       fC[7][0] = 0.005;
4386       fC[7][1] = 0.005;
4387       fC[7][2] = 1.0000;
4388       fC[7][3] = 0.1519;
4389       fC[7][4] = 0.1989;
4390
4391       fC[8][0] = 0.005;
4392       fC[8][1] = 0.005;
4393       fC[8][2] = 1.0000;
4394       fC[8][3] = 0.1817;
4395       fC[8][4] = 0.2472;
4396
4397       fC[9][0] = 0.005;
4398       fC[9][1] = 0.005;
4399       fC[9][2] = 1.0000;
4400       fC[9][3] = 0.2429;
4401       fC[9][4] = 0.2684;
4402
4403       fC[10][0] = 0.005;
4404       fC[10][1] = 0.005;
4405       fC[10][2] = 1.0000;
4406       fC[10][3] = 0.2760;
4407       fC[10][4] = 0.3098;
4408
4409       fC[11][0] = 0.005;
4410       fC[11][1] = 0.005;
4411       fC[11][2] = 1.0000;
4412       fC[11][3] = 0.2673;
4413       fC[11][4] = 0.3198;
4414
4415       fC[12][0] = 0.005;
4416       fC[12][1] = 0.005;
4417       fC[12][2] = 1.0000;
4418       fC[12][3] = 0.3165;
4419       fC[12][4] = 0.3564;
4420
4421       fC[13][0] = 0.005;
4422       fC[13][1] = 0.005;
4423       fC[13][2] = 1.0000;
4424       fC[13][3] = 0.3526;
4425       fC[13][4] = 0.3011;
4426
4427       fC[14][0] = 0.005;
4428       fC[14][1] = 0.005;
4429       fC[14][2] = 1.0000;
4430       fC[14][3] = 0.3788;
4431       fC[14][4] = 0.3011;
4432
4433       fC[15][0] = 0.005;
4434       fC[15][1] = 0.005;
4435       fC[15][2] = 1.0000;
4436       fC[15][3] = 0.3788;
4437       fC[15][4] = 0.3011;
4438
4439       fC[16][0] = 0.005;
4440       fC[16][1] = 0.005;
4441       fC[16][2] = 1.0000;
4442       fC[16][3] = 0.3788;
4443       fC[16][4] = 0.3011;
4444
4445       fC[17][0] = 0.005;
4446       fC[17][1] = 0.005;
4447       fC[17][2] = 1.0000;
4448       fC[17][3] = 0.3788;
4449       fC[17][4] = 0.3011;
4450   }
4451
4452   for(Int_t i=18;i<fgkPIDptBin;i++){
4453     fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
4454     fC[i][0] = fC[17][0];
4455     fC[i][1] = fC[17][1];
4456     fC[i][2] = fC[17][2];
4457     fC[i][3] = fC[17][3];
4458     fC[i][4] = fC[17][4];
4459   }
4460 }
4461
4462 //---------------------------------------------------------------//
4463 Bool_t AliFlowTrackCuts::TPCTOFagree(const AliVTrack *track)
4464 {
4465   //check pid agreement between TPC and TOF
4466   Bool_t status = kFALSE;
4467
4468   const Float_t c = 2.99792457999999984e-02;
4469
4470   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
4471
4472
4473   Double_t exptimes[9];
4474   track->GetIntegratedTimes(exptimes);
4475
4476   Float_t dedx = track->GetTPCsignal();
4477
4478   Float_t p = track->P();
4479   Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
4480   Float_t tl = exptimes[0]*c; // to work both for ESD and AOD
4481
4482   Float_t betagammares =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4483
4484   Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
4485
4486 //  printf("betagamma1 = %f\n",betagamma1);
4487
4488   if(betagamma1 < 0.1) betagamma1 = 0.1;
4489
4490   if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
4491   else betagamma1 = 100;
4492
4493   Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
4494 //  printf("betagamma2 = %f\n",betagamma2);
4495
4496   if(betagamma2 < 0.1) betagamma2 = 0.1;
4497
4498   if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
4499   else betagamma2 = 100;
4500
4501
4502   Float_t momtpc=track->GetTPCmomentum();
4503
4504   for(Int_t i=0;i < 5;i++){
4505     Float_t resolutionTOF =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
4506     if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
4507       Float_t dedxExp = 0;
4508       if(i==0) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
4509       else if(i==1) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
4510       else if(i==2) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
4511       else if(i==3) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
4512       else if(i==4) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
4513
4514       Float_t resolutionTPC = 2;
4515       if(i==0) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
4516       else if(i==1) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
4517       else if(i==2) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
4518       else if(i==3) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
4519       else if(i==4) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4520
4521       if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
4522         status = kTRUE;
4523       }
4524     }
4525   }
4526
4527   Float_t bb1 =  fESDpid.GetTPCResponse().Bethe(betagamma1);
4528   Float_t bb2 =  fESDpid.GetTPCResponse().Bethe(betagamma2);
4529   Float_t bbM =  fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
4530
4531
4532   //  status = kFALSE;
4533   // for nuclei
4534   Float_t resolutionTOFpr =   fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4535   Float_t resolutionTPCpr =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4536   if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4537      status = kTRUE;
4538   }
4539   else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4540      status = kTRUE;
4541   }
4542   else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4543      status = kTRUE;
4544   }
4545
4546   return status;
4547 }
4548 // end part added by F. Noferini
4549 //-----------------------------------------------------------------------
4550
4551 //-----------------------------------------------------------------------
4552 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
4553 {
4554   //get the name of the particle id source
4555   switch (s)
4556   {
4557     case kTPCdedx:
4558       return "TPCdedx";
4559     case kTPCbayesian:
4560       return "TPCbayesian";
4561     case kTOFbeta:
4562       return "TOFbeta";
4563     case kTPCpid:
4564       return "TPCpid";
4565     case kTOFpid:
4566       return "TOFpid";
4567     case kTOFbayesian:
4568       return "TOFbayesianPID";
4569     case kTOFbetaSimple:
4570       return "TOFbetaSimple";
4571     case kTPCNuclei:
4572       return "TPCnuclei";
4573     default:
4574       return "NOPID";
4575   }
4576 }
4577
4578 //-----------------------------------------------------------------------
4579 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
4580 {
4581   //return the name of the selected parameter type
4582   switch (type)
4583   {
4584     case kMC:
4585       return "MC";
4586     case kGlobal:
4587       return "Global";
4588     case kTPCstandalone:
4589       return "TPCstandalone";
4590     case kSPDtracklet:
4591       return "SPDtracklets";
4592     case kPMD:
4593       return "PMD";
4594     case kVZERO:
4595       return "VZERO";
4596     case kMUON:       // XZhang 20120604
4597       return "MUON";  // XZhang 20120604
4598     case kKink:
4599       return "Kink";
4600     case kV0:
4601       return "V0";
4602     default:
4603       return "unknown";
4604   }
4605 }
4606
4607 //-----------------------------------------------------------------------
4608 Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
4609 {
4610   //check PMD specific cuts
4611   //clean up from last iteration, and init label
4612   Int_t   det   = track->GetDetector();
4613   //Int_t   smn   = track->GetSmn();
4614   Float_t clsX  = track->GetClusterX();
4615   Float_t clsY  = track->GetClusterY();
4616   Float_t clsZ  = track->GetClusterZ();
4617   Float_t ncell = track->GetClusterCells();
4618   Float_t adc   = track->GetClusterADC();
4619
4620   fTrack = NULL;
4621   fMCparticle=NULL;
4622   fTrackLabel=-996;
4623
4624   fTrackEta = GetPmdEta(clsX,clsY,clsZ);
4625   fTrackPhi = GetPmdPhi(clsX,clsY);
4626   fTrackWeight = 1.0;
4627
4628   Bool_t pass=kTRUE;
4629   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
4630   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
4631   if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
4632   if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
4633   if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
4634
4635   return pass;
4636 }
4637
4638 //-----------------------------------------------------------------------
4639 Bool_t AliFlowTrackCuts::PassesVZEROcuts(Int_t id)
4640 {
4641   //check VZERO cuts
4642   if (id<0) return kFALSE;
4643
4644   //clean up from last iter
4645   ClearTrack();
4646
4647   fTrackPhi = TMath::PiOver4()*(0.5+id%8);
4648
4649   // 10102013 weighting vzero tiles - rbertens@cern.ch
4650   if(!fVZEROgainEqualization) {
4651       // if for some reason the equalization is not initialized (e.g. 2011 data)
4652       // the fVZEROxpol[] weights are used to enable or disable vzero rings
4653     if(id<32) {   // v0c side
4654       fTrackEta = -3.45+0.5*(id/8);
4655       if(id < 8) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[0];
4656       else if (id < 16 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[1];
4657       else if (id < 24 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[2];
4658       else if (id < 32 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[3];
4659     } else {       // v0a side
4660       fTrackEta = +4.8-0.6*((id/8)-4);
4661       if( id < 40) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[0];
4662       else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[1];
4663       else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[2];
4664       else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[3];
4665     }
4666   } else { // the equalization is initialized
4667      // note that disabled rings have already been excluded on calibration level in
4668      // AliFlowEvent (so for a disabled ring, fVZEROxpol is zero
4669     if(id<32) {    // v0c side
4670       fTrackEta = -3.45+0.5*(id/8);
4671       if(id < 8) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4672       else if (id < 16 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4673       else if (id < 24 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4674       else if (id < 32 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[3]/fVZEROgainEqualization->GetBinContent(1+id);
4675     } else {       // v0a side
4676       fTrackEta = +4.8-0.6*((id/8)-4);
4677       if( id < 40) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4678       else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4679       else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4680       else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[3]/fVZEROgainEqualization->GetBinContent(1+id);
4681     }
4682     // printf ( " tile %i and weight %.2f \n", id, fTrackWeight);
4683   }
4684
4685   if (fLinearizeVZEROresponse && id < 64)
4686   {
4687     //this is only needed in pass1 of LHC10h
4688     Float_t multVZERO[fgkNumberOfVZEROtracks];
4689     Float_t dummy=0.0;
4690     AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multVZERO);
4691     fTrackWeight = multVZERO[id];
4692   }
4693
4694   Bool_t pass=kTRUE;
4695   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
4696   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
4697
4698   return pass;
4699 }
4700
4701 //-----------------------------------------------------------------------------
4702 Bool_t AliFlowTrackCuts::PassesMuonCuts(AliVParticle* vparticle)
4703 {
4704 // XZhang 20120604
4705   ClearTrack();
4706   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
4707   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
4708   else fMCparticle=NULL;
4709
4710   AliESDMuonTrack *esdTrack = dynamic_cast<AliESDMuonTrack*>(vparticle);
4711   AliAODTrack     *aodTrack = dynamic_cast<AliAODTrack*>(vparticle);
4712   if ((!esdTrack) && (!aodTrack)) return kFALSE;
4713   if (!fMuonTrackCuts->IsSelected(vparticle)) return kFALSE;
4714   HandleVParticle(vparticle);   if (!fTrack)  return kFALSE;
4715   return kTRUE;
4716 }
4717
4718 //----------------------------------------------------------------------------//
4719 Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
4720 {
4721   //get PMD "track" eta coordinate
4722   Float_t rpxpy, theta, eta;
4723   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
4724   theta  = TMath::ATan2(rpxpy,zPos);
4725   eta    = -TMath::Log(TMath::Tan(0.5*theta));
4726   return eta;
4727 }
4728
4729 //--------------------------------------------------------------------------//
4730 Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
4731 {
4732   //Get PMD "track" phi coordinate
4733   Float_t pybypx, phi = 0., phi1;
4734   if(xPos==0)
4735     {
4736       if(yPos>0) phi = 90.;
4737       if(yPos<0) phi = 270.;
4738     }
4739   if(xPos != 0)
4740     {
4741       pybypx = yPos/xPos;
4742       if(pybypx < 0) pybypx = - pybypx;
4743       phi1 = TMath::ATan(pybypx)*180./3.14159;
4744
4745       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
4746       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
4747       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
4748       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
4749
4750     }
4751   phi = phi*3.14159/180.;
4752   return   phi;
4753 }
4754
4755 //---------------------------------------------------------------//
4756 void AliFlowTrackCuts::Browse(TBrowser* b)
4757 {
4758   //some browsing capabilities
4759   if (fQA) b->Add(fQA);
4760 }
4761
4762 //---------------------------------------------------------------//
4763 Long64_t AliFlowTrackCuts::Merge(TCollection* list)
4764 {
4765   //merge
4766   if (!fQA || !list) return 0;
4767   if (list->IsEmpty()) return 0;
4768   AliFlowTrackCuts* obj=NULL;
4769   TList tmplist;
4770   TIter next(list);
4771   while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
4772   {
4773     if (obj==this) continue;
4774     tmplist.Add(obj->GetQA());
4775   }
4776   return fQA->Merge(&tmplist);
4777 }
4778
4779
4780