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   AliVVZERO* vvzero = dynamic_cast<AliVVZERO*>(obj);    // downcast to base class
 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 (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1182     {
1183       if (!PassesAODpidCut(track)) pass=kFALSE;
1184     }
1185
1186   if (fQA) {
1187     // changed 04062014 used to be filled before possible PID cut
1188     Double_t momTPC = track->GetTPCmomentum();
1189     QAbefore( 1)->Fill(momTPC,dedx);
1190     QAbefore( 5)->Fill(track->Pt(),track->DCA());
1191     QAbefore( 6)->Fill(track->Pt(),track->ZAtDCA());
1192     if (pass) QAafter( 1)->Fill(momTPC,dedx);
1193     if (pass) QAafter( 5)->Fill(track->Pt(),track->DCA());
1194     if (pass) QAafter( 6)->Fill(track->Pt(),track->ZAtDCA());
1195     QAbefore( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1196     if (pass) QAafter(  8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1197     QAbefore( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1198     if (pass) QAafter(  9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1199     QAbefore(10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1200     if (pass) QAafter( 10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1201     QAbefore(11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1202     if (pass) QAafter( 11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1203     QAbefore(12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1204     if (pass) QAafter( 12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1205   }
1206
1207
1208   return pass;
1209 }
1210
1211 //_______________________________________________________________________
1212 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
1213 {
1214   //check cuts on ESD tracks
1215   Bool_t pass=kTRUE;
1216   Float_t dcaxy=0.0;
1217   Float_t dcaz=0.0;
1218   track->GetImpactParameters(dcaxy,dcaz);
1219   const AliExternalTrackParam* pout = track->GetOuterParam();
1220   const AliExternalTrackParam* pin = track->GetInnerParam();
1221   if (fIgnoreTPCzRange)
1222   {
1223     if (pin&&pout)
1224     {
1225       Double_t zin = pin->GetZ();
1226       Double_t zout = pout->GetZ();
1227       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
1228       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
1229       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
1230     }
1231   }
1232
1233   Int_t ntpccls = ( fParamType==kTPCstandalone )?
1234                     track->GetTPCNclsIter1():track->GetTPCNcls();
1235   if (fCutChi2PerClusterTPC)
1236   {
1237     Float_t tpcchi2 = (fParamType==kTPCstandalone)?
1238                        track->GetTPCchi2Iter1():track->GetTPCchi2();
1239     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
1240     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
1241       pass=kFALSE;
1242   }
1243
1244   if (fCutMinimalTPCdedx)
1245   {
1246     if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
1247   }
1248
1249   if (fCutNClustersTPC)
1250   {
1251     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1252   }
1253
1254   Int_t nitscls = track->GetNcls(0);
1255   if (fCutNClustersITS)
1256   {
1257     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1258   }
1259
1260   //some stuff is still handled by AliESDtrackCuts class - delegate
1261   if (fAliESDtrackCuts)
1262   {
1263     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
1264   }
1265
1266   //PID part with pid QA
1267   Double_t beta = GetBeta(track);
1268   Double_t dedx = Getdedx(track);
1269   if (fQA)
1270   {
1271     if (pass) QAbefore(0)->Fill(track->GetP(),beta);
1272     if (pass) QAbefore(1)->Fill(pin->GetP(),dedx);
1273   }
1274   if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1275   {
1276     if (!PassesESDpidCut(track)) pass=kFALSE;
1277   }
1278   if (fCutRejectElectronsWithTPCpid)
1279   {
1280     //reject electrons using standard TPC pid
1281     //TODO this should be rethought....
1282     Double_t pidTPC[AliPID::kSPECIES];
1283     track->GetTPCpid(pidTPC);
1284     if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
1285   }
1286   if (fQA)
1287   {
1288     if (pass) QAafter(0)->Fill(track->GetP(),beta);
1289     if (pass) QAafter(1)->Fill(pin->GetP(),dedx);
1290   }
1291   //end pid part with qa
1292
1293   if (fQA)
1294   {
1295     Double_t pt = track->Pt();
1296     QAbefore(5)->Fill(pt,dcaxy);
1297     QAbefore(6)->Fill(pt,dcaz);
1298     if (pass) QAafter(5)->Fill(pt,dcaxy);
1299     if (pass) QAafter(6)->Fill(pt,dcaz);
1300     QAbefore(17)->Fill(Float_t(track->GetKinkIndex(0)));
1301     if (pass) QAafter(17)->Fill(Float_t(track->GetKinkIndex(0)));
1302   }
1303
1304   return pass;
1305 }
1306
1307 //-----------------------------------------------------------------------
1308 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
1309 {
1310   //handle the general case
1311   switch (fParamType)
1312   {
1313     default:
1314       fTrack = track;
1315       break;
1316   }
1317 }
1318
1319 //-----------------------------------------------------------------------
1320 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
1321 {
1322   //handle esd track
1323   switch (fParamType)
1324   {
1325     case kGlobal:
1326       fTrack = track;
1327       break;
1328     case kTPCstandalone:
1329       if (!track->FillTPCOnlyTrack(fTPCtrack))
1330       {
1331         fTrack=NULL;
1332         fMCparticle=NULL;
1333         fTrackLabel=-998;
1334         return;
1335       }
1336       fTrack = &fTPCtrack;
1337       //recalculate the label and mc particle, they may differ as TPClabel != global label
1338       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1339       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1340       else fMCparticle=NULL;
1341       break;
1342     default:
1343       if (fForceTPCstandalone)
1344       {
1345         if (!track->FillTPCOnlyTrack(fTPCtrack))
1346         {
1347           fTrack=NULL;
1348           fMCparticle=NULL;
1349           fTrackLabel=-998;
1350           return;
1351         }
1352         fTrack = &fTPCtrack;
1353         //recalculate the label and mc particle, they may differ as TPClabel != global label
1354         fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1355         if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1356         else fMCparticle=NULL;
1357       }
1358       else
1359         fTrack = track;
1360       break;
1361   }
1362 }
1363
1364 //-----------------------------------------------------------------------
1365 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
1366 {
1367   //calculate the number of track in given event.
1368   //if argument is NULL(default) take the event attached
1369   //by SetEvent()
1370   Int_t multiplicity = 0;
1371   if (!event)
1372   {
1373     for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
1374     {
1375       if (IsSelected(GetInputObject(i))) multiplicity++;
1376     }
1377   }
1378   else
1379   {
1380     for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
1381     {
1382       if (IsSelected(event->GetTrack(i))) multiplicity++;
1383     }
1384   }
1385   return multiplicity;
1386 }
1387
1388 //-----------------------------------------------------------------------
1389 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts()
1390 {
1391   //returns the lhc10h vzero track cuts, this function
1392   //is left here for backward compatibility
1393   //if a run is recognized as 11h, the calibration method will
1394   //switch to 11h calbiration, which means that the cut
1395   //object is updated but not replaced.
1396   //calibratin is only available for PbPb runs
1397   return GetStandardVZEROOnlyTrackCuts2010();
1398 }
1399 //-----------------------------------------------------------------------
1400 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010()
1401 {
1402   //get standard VZERO cuts
1403   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2010");
1404   cuts->SetParamType(kVZERO);
1405   cuts->SetEtaRange( -10, +10 );
1406   cuts->SetPhiMin( 0 );
1407   cuts->SetPhiMax( TMath::TwoPi() );
1408   // options for the reweighting
1409   cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1410   cuts->SetApplyRecentering(kTRUE);
1411   // to exclude a ring , do e.g.
1412   // cuts->SetUseVZERORing(7, kFALSE);
1413   return cuts;
1414 }
1415 //-----------------------------------------------------------------------
1416 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011()
1417 {
1418   //get standard VZERO cuts for 2011 data
1419   //in this case, the vzero segments will be weighted by
1420   //VZEROEqMultiplicity,
1421   //if recentering is enableded, the sub-q vectors
1422   //will be taken from the event header, so make sure to run
1423   //the VZERO event plane selection task before this task !
1424   //recentering replaces the already evaluated q-vectors, so
1425   //when chosen, additional settings (e.g. excluding rings)
1426   //have no effect. recentering is true by default
1427   //
1428   //NOTE user is responsible for running the vzero event plane
1429   //selection task in advance, e.g. add to your launcher macro
1430   //
1431   //  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C");
1432   //  AddTaskVZEROEPSelection();
1433   //
1434   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2011");
1435   cuts->SetParamType(kVZERO);
1436   cuts->SetEtaRange( -10, +10 );
1437   cuts->SetPhiMin( 0 );
1438   cuts->SetPhiMax( TMath::TwoPi() );
1439   cuts->SetApplyRecentering(kTRUE);
1440   cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1441  return cuts;
1442 }
1443 //-----------------------------------------------------------------------
1444 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
1445 {
1446   //get standard cuts
1447   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global tracks");
1448   cuts->SetParamType(kGlobal);
1449   cuts->SetPtRange(0.2,5.);
1450   cuts->SetEtaRange(-0.8,0.8);
1451   cuts->SetMinNClustersTPC(70);
1452   cuts->SetMinChi2PerClusterTPC(0.1);
1453   cuts->SetMaxChi2PerClusterTPC(4.0);
1454   cuts->SetMinNClustersITS(2);
1455   cuts->SetRequireITSRefit(kTRUE);
1456   cuts->SetRequireTPCRefit(kTRUE);
1457   cuts->SetMaxDCAToVertexXY(0.3);
1458   cuts->SetMaxDCAToVertexZ(0.3);
1459   cuts->SetAcceptKinkDaughters(kFALSE);
1460   cuts->SetMinimalTPCdedx(10.);
1461
1462   return cuts;
1463 }
1464
1465 //-----------------------------------------------------------------------
1466 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()
1467 {
1468   //get standard cuts
1469   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone 2010");
1470   cuts->SetParamType(kTPCstandalone);
1471   cuts->SetPtRange(0.2,5.);
1472   cuts->SetEtaRange(-0.8,0.8);
1473   cuts->SetMinNClustersTPC(70);
1474   cuts->SetMinChi2PerClusterTPC(0.2);
1475   cuts->SetMaxChi2PerClusterTPC(4.0);
1476   cuts->SetMaxDCAToVertexXY(3.0);
1477   cuts->SetMaxDCAToVertexZ(3.0);
1478   cuts->SetDCAToVertex2D(kTRUE);
1479   cuts->SetAcceptKinkDaughters(kFALSE);
1480   cuts->SetMinimalTPCdedx(10.);
1481
1482   return cuts;
1483 }
1484
1485 //-----------------------------------------------------------------------
1486 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts()
1487 {
1488   //get standard cuts
1489   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone");
1490   cuts->SetParamType(kTPCstandalone);
1491   cuts->SetPtRange(0.2,5.);
1492   cuts->SetEtaRange(-0.8,0.8);
1493   cuts->SetMinNClustersTPC(70);
1494   cuts->SetMinChi2PerClusterTPC(0.2);
1495   cuts->SetMaxChi2PerClusterTPC(4.0);
1496   cuts->SetMaxDCAToVertexXY(3.0);
1497   cuts->SetMaxDCAToVertexZ(3.0);
1498   cuts->SetDCAToVertex2D(kTRUE);
1499   cuts->SetAcceptKinkDaughters(kFALSE);
1500   cuts->SetMinimalTPCdedx(10.);
1501
1502   return cuts;
1503 }
1504
1505 //-----------------------------------------------------------------------
1506 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
1507 {
1508   //get standard cuts
1509   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
1510   delete cuts->fAliESDtrackCuts;
1511   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
1512   cuts->SetParamType(kGlobal);
1513   return cuts;
1514 }
1515
1516 //-----------------------------------------------------------------------------
1517 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardMuonTrackCuts(Bool_t isMC, Int_t passN)
1518 {
1519 // XZhang 20120604
1520   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard muon track cuts");
1521   cuts->SetParamType(kMUON);
1522   cuts->SetStandardMuonTrackCuts();
1523   cuts->SetIsMuonMC(isMC);
1524   cuts->SetMuonPassNumber(passN);
1525   return cuts;
1526 }
1527
1528 //-----------------------------------------------------------------------
1529 //AliFlowTrack* AliFlowTrackCuts::FillFlowTrackV0(TObjArray* trackCollection, Int_t trackIndex) const
1530 //{
1531 //  //fill flow track from a reconstructed V0 (topological)
1532 //  AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1533 //  if (flowtrack)
1534 //  {
1535 //    flowtrack->Clear();
1536 //  }
1537 //  else
1538 //  {
1539 //    flowtrack = new AliFlowCandidateTrack();
1540 //    trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1541 //  }
1542 //
1543 //  TParticle *tmpTParticle=NULL;
1544 //  AliMCParticle* tmpAliMCParticle=NULL;
1545 //  AliExternalTrackParam* externalParams=NULL;
1546 //  AliESDtrack* esdtrack=NULL;
1547 //  switch(fParamMix)
1548 //  {
1549 //    case kPure:
1550 //      flowtrack->Set(fTrack);
1551 //      break;
1552 //    case kTrackWithMCkine:
1553 //      flowtrack->Set(fMCparticle);
1554 //      break;
1555 //    case kTrackWithMCPID:
1556 //      flowtrack->Set(fTrack);
1557 //      //flowtrack->setPID(...) from mc, when implemented
1558 //      break;
1559 //    case kTrackWithMCpt:
1560 //      if (!fMCparticle) return NULL;
1561 //      flowtrack->Set(fTrack);
1562 //      flowtrack->SetPt(fMCparticle->Pt());
1563 //      break;
1564 //    case kTrackWithPtFromFirstMother:
1565 //      if (!fMCparticle) return NULL;
1566 //      flowtrack->Set(fTrack);
1567 //      tmpTParticle = fMCparticle->Particle();
1568 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1569 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
1570 //      break;
1571 //    case kTrackWithTPCInnerParams:
1572 //      esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1573 //      if (!esdtrack) return NULL;
1574 //      externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1575 //      if (!externalParams) return NULL;
1576 //      flowtrack->Set(externalParams);
1577 //      break;
1578 //    default:
1579 //      flowtrack->Set(fTrack);
1580 //      break;
1581 //  }
1582 //  if (fParamType==kMC)
1583 //  {
1584 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
1585 //    flowtrack->SetID(fTrackLabel);
1586 //  }
1587 //  else if (dynamic_cast<AliAODTrack*>(fTrack))
1588 //  {
1589 //    if (fParamType==kMUON)                            // XZhang 20120604
1590 //      flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1591 //    else                                              // XZhang 20120604
1592 //      flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1593 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1594 //  }
1595 //  else if (dynamic_cast<AliMCParticle*>(fTrack))
1596 //  {
1597 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
1598 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1599 //  }
1600 //
1601 //  if (fV0)
1602 //  {
1603 //    //add daughter indices
1604 //  }
1605 //
1606 //  return flowtrack;
1607 //}
1608
1609 //-----------------------------------------------------------------------
1610 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackKink(TObjArray* trackCollection, Int_t trackIndex) const
1611 {
1612   //fill flow track from AliVParticle (ESD,AOD,MC)
1613   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1614   if (flowtrack)
1615   {
1616     flowtrack->Clear();
1617   }
1618   else
1619   {
1620     flowtrack = new AliFlowCandidateTrack();
1621     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1622   }
1623
1624   TParticle *tmpTParticle=NULL;
1625   AliMCParticle* tmpAliMCParticle=NULL;
1626   AliExternalTrackParam* externalParams=NULL;
1627   AliESDtrack* esdtrack=NULL;
1628   switch(fParamMix)
1629   {
1630     case kPure:
1631       flowtrack->Set(fTrack);
1632       break;
1633     case kTrackWithMCkine:
1634       flowtrack->Set(fMCparticle);
1635       break;
1636     case kTrackWithMCPID:
1637       flowtrack->Set(fTrack);
1638       //flowtrack->setPID(...) from mc, when implemented
1639       break;
1640     case kTrackWithMCpt:
1641       if (!fMCparticle) return NULL;
1642       flowtrack->Set(fTrack);
1643       flowtrack->SetPt(fMCparticle->Pt());
1644       break;
1645     case kTrackWithPtFromFirstMother:
1646       if (!fMCparticle) return NULL;
1647       flowtrack->Set(fTrack);
1648       tmpTParticle = fMCparticle->Particle();
1649       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1650       flowtrack->SetPt(tmpAliMCParticle->Pt());
1651       break;
1652     case kTrackWithTPCInnerParams:
1653       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1654       if (!esdtrack) return NULL;
1655       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1656       if (!externalParams) return NULL;
1657       flowtrack->Set(externalParams);
1658       break;
1659     default:
1660       flowtrack->Set(fTrack);
1661       break;
1662   }
1663   if (fParamType==kMC)
1664   {
1665     flowtrack->SetSource(AliFlowTrack::kFromMC);
1666     flowtrack->SetID(fTrackLabel);
1667   }
1668   else if (dynamic_cast<AliESDtrack*>(fTrack))
1669   {
1670     flowtrack->SetSource(AliFlowTrack::kFromESD);
1671     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1672   }
1673   else if (dynamic_cast<AliESDMuonTrack*>(fTrack))                                  // XZhang 20120604
1674   {                                                                                 // XZhang 20120604
1675     flowtrack->SetSource(AliFlowTrack::kFromMUON);                                  // XZhang 20120604
1676     flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID());  // XZhang 20120604
1677   }                                                                                 // XZhang 20120604
1678   else if (dynamic_cast<AliAODTrack*>(fTrack))
1679   {
1680     if (fParamType==kMUON)                            // XZhang 20120604
1681       flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1682     else                                              // XZhang 20120604
1683       flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1684     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1685   }
1686   else if (dynamic_cast<AliMCParticle*>(fTrack))
1687   {
1688     flowtrack->SetSource(AliFlowTrack::kFromMC);
1689     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1690   }
1691
1692   if (fKink)
1693   {
1694     Int_t indexMother = fKink->GetIndex(0);
1695     Int_t indexDaughter = fKink->GetIndex(1);
1696     flowtrack->SetID(indexMother);
1697     flowtrack->AddDaughter(indexDaughter);
1698     flowtrack->SetMass(1.);
1699     flowtrack->SetSource(AliFlowTrack::kFromKink);
1700   }
1701
1702   flowtrack->SetMass(fTrackMass);
1703
1704   return flowtrack;
1705 }
1706
1707 //-----------------------------------------------------------------------
1708 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackGeneric(TObjArray* trackCollection, Int_t trackIndex) const
1709 {
1710   //fill a flow track from tracklet,vzero,pmd,...
1711   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1712   if (flowtrack)
1713   {
1714     flowtrack->Clear();
1715   }
1716   else
1717   {
1718     flowtrack = new AliFlowTrack();
1719     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1720   }
1721
1722   if (FillFlowTrackGeneric(flowtrack)) return flowtrack;
1723   else
1724   {
1725     trackCollection->RemoveAt(trackIndex);
1726     delete flowtrack;
1727     return NULL;
1728   }
1729 }
1730
1731 //-----------------------------------------------------------------------
1732 Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
1733 {
1734   //fill a flow track from tracklet,vzero,pmd,...
1735   TParticle *tmpTParticle=NULL;
1736   AliMCParticle* tmpAliMCParticle=NULL;
1737   switch (fParamMix)
1738   {
1739     case kPure:
1740       flowtrack->SetPhi(fTrackPhi);
1741       flowtrack->SetEta(fTrackEta);
1742       flowtrack->SetWeight(fTrackWeight);
1743       flowtrack->SetPt(fTrackPt);
1744       flowtrack->SetMass(fTrackMass);
1745       break;
1746     case kTrackWithMCkine:
1747       if (!fMCparticle) return kFALSE;
1748       flowtrack->SetPhi( fMCparticle->Phi() );
1749       flowtrack->SetEta( fMCparticle->Eta() );
1750       flowtrack->SetPt( fMCparticle->Pt() );
1751       flowtrack->SetMass(fTrackMass);
1752       break;
1753     case kTrackWithMCpt:
1754       if (!fMCparticle) return kFALSE;
1755       flowtrack->SetPhi(fTrackPhi);
1756       flowtrack->SetEta(fTrackEta);
1757       flowtrack->SetWeight(fTrackWeight);
1758       flowtrack->SetPt(fMCparticle->Pt());
1759       flowtrack->SetMass(fTrackMass);
1760       break;
1761     case kTrackWithPtFromFirstMother:
1762       if (!fMCparticle) return kFALSE;
1763       flowtrack->SetPhi(fTrackPhi);
1764       flowtrack->SetEta(fTrackEta);
1765       flowtrack->SetWeight(fTrackWeight);
1766       tmpTParticle = fMCparticle->Particle();
1767       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1768       flowtrack->SetPt(tmpAliMCParticle->Pt());
1769       flowtrack->SetMass(fTrackMass);
1770       break;
1771     default:
1772       flowtrack->SetPhi(fTrackPhi);
1773       flowtrack->SetEta(fTrackEta);
1774       flowtrack->SetWeight(fTrackWeight);
1775       flowtrack->SetPt(fTrackPt);
1776       flowtrack->SetMass(fTrackMass);
1777       break;
1778   }
1779   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1780   return kTRUE;
1781 }
1782
1783 //-----------------------------------------------------------------------
1784 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackVParticle(TObjArray* trackCollection, Int_t trackIndex) const
1785 {
1786   //fill flow track from AliVParticle (ESD,AOD,MC)
1787
1788   AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1789   if (flowtrack)
1790   {
1791     flowtrack->Clear();
1792   }
1793   else
1794   {
1795     flowtrack = new AliFlowTrack();
1796     trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1797   }
1798
1799   if (FillFlowTrackVParticle(flowtrack)) return flowtrack;
1800   else
1801   {
1802     trackCollection->RemoveAt(trackIndex);
1803     delete flowtrack;
1804     return NULL;
1805   }
1806 }
1807
1808 //-----------------------------------------------------------------------
1809 Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
1810 {
1811   //fill flow track from AliVParticle (ESD,AOD,MC)
1812   if (!fTrack) return kFALSE;
1813   TParticle *tmpTParticle=NULL;
1814   AliMCParticle* tmpAliMCParticle=NULL;
1815   AliExternalTrackParam* externalParams=NULL;
1816   AliESDtrack* esdtrack=NULL;
1817   switch(fParamMix)
1818   {
1819     case kPure:
1820       flowtrack->Set(fTrack);
1821       break;
1822     case kTrackWithMCkine:
1823       flowtrack->Set(fMCparticle);
1824       break;
1825     case kTrackWithMCPID:
1826       flowtrack->Set(fTrack);
1827       //flowtrack->setPID(...) from mc, when implemented
1828       break;
1829     case kTrackWithMCpt:
1830       if (!fMCparticle) return kFALSE;
1831       flowtrack->Set(fTrack);
1832       flowtrack->SetPt(fMCparticle->Pt());
1833       break;
1834     case kTrackWithPtFromFirstMother:
1835       if (!fMCparticle) return kFALSE;
1836       flowtrack->Set(fTrack);
1837       tmpTParticle = fMCparticle->Particle();
1838       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1839       flowtrack->SetPt(tmpAliMCParticle->Pt());
1840       break;
1841     case kTrackWithTPCInnerParams:
1842       esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1843       if (!esdtrack) return kFALSE;
1844       externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1845       if (!externalParams) return kFALSE;
1846       flowtrack->Set(externalParams);
1847       break;
1848     default:
1849       flowtrack->Set(fTrack);
1850       break;
1851   }
1852   if (fParamType==kMC)
1853   {
1854     flowtrack->SetSource(AliFlowTrack::kFromMC);
1855     flowtrack->SetID(fTrackLabel);
1856   }
1857   else if (dynamic_cast<AliESDtrack*>(fTrack))
1858   {
1859     flowtrack->SetSource(AliFlowTrack::kFromESD);
1860     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1861   }
1862   else if (dynamic_cast<AliESDMuonTrack*>(fTrack))                                  // XZhang 20120604
1863   {                                                                                 // XZhang 20120604
1864     flowtrack->SetSource(AliFlowTrack::kFromMUON);                                  // XZhang 20120604
1865     flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID());  // XZhang 20120604
1866   }                                                                                 // XZhang 20120604
1867   else if (dynamic_cast<AliAODTrack*>(fTrack))
1868   {
1869     if (fParamType==kMUON)                            // XZhang 20120604
1870       flowtrack->SetSource(AliFlowTrack::kFromMUON);  // XZhang 20120604
1871     else                                              // XZhang 20120604
1872       flowtrack->SetSource(AliFlowTrack::kFromAOD);   // XZhang 20120604
1873     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1874   }
1875   else if (dynamic_cast<AliMCParticle*>(fTrack))
1876   {
1877     flowtrack->SetSource(AliFlowTrack::kFromMC);
1878     flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1879   }
1880   flowtrack->SetMass(fTrackMass);
1881   return kTRUE;
1882 }
1883
1884 //-----------------------------------------------------------------------
1885 AliFlowTrack* AliFlowTrackCuts::FillFlowTrack(TObjArray* trackCollection, Int_t trackIndex) const
1886 {
1887   //fill a flow track constructed from whatever we applied cuts on
1888   //return true on success
1889   switch (fParamType)
1890   {
1891     case kSPDtracklet:
1892       return FillFlowTrackGeneric(trackCollection, trackIndex);
1893     case kPMD:
1894       return FillFlowTrackGeneric(trackCollection, trackIndex);
1895     case kVZERO:
1896       return FillFlowTrackGeneric(trackCollection, trackIndex);
1897     case kKink:
1898       return FillFlowTrackKink(trackCollection, trackIndex);
1899     //case kV0:
1900     //  return FillFlowTrackV0(trackCollection, trackIndex);
1901     default:
1902       return FillFlowTrackVParticle(trackCollection, trackIndex);
1903   }
1904 }
1905
1906 //-----------------------------------------------------------------------
1907 Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
1908 {
1909   //fill a flow track constructed from whatever we applied cuts on
1910   //return true on success
1911   switch (fParamType)
1912   {
1913     case kSPDtracklet:
1914       return FillFlowTrackGeneric(track);
1915     case kPMD:
1916       return FillFlowTrackGeneric(track);
1917     case kVZERO:
1918       return FillFlowTrackGeneric(track);
1919     default:
1920       return FillFlowTrackVParticle(track);
1921   }
1922 }
1923
1924 ////-----------------------------------------------------------------------
1925 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
1926 //{
1927 //  //make a flow track from tracklet
1928 //  AliFlowTrack* flowtrack=NULL;
1929 //  TParticle *tmpTParticle=NULL;
1930 //  AliMCParticle* tmpAliMCParticle=NULL;
1931 //  switch (fParamMix)
1932 //  {
1933 //    case kPure:
1934 //      flowtrack = new AliFlowTrack();
1935 //      flowtrack->SetPhi(fTrackPhi);
1936 //      flowtrack->SetEta(fTrackEta);
1937 //      flowtrack->SetWeight(fTrackWeight);
1938 //      break;
1939 //    case kTrackWithMCkine:
1940 //      if (!fMCparticle) return NULL;
1941 //      flowtrack = new AliFlowTrack();
1942 //      flowtrack->SetPhi( fMCparticle->Phi() );
1943 //      flowtrack->SetEta( fMCparticle->Eta() );
1944 //      flowtrack->SetPt( fMCparticle->Pt() );
1945 //      break;
1946 //    case kTrackWithMCpt:
1947 //      if (!fMCparticle) return NULL;
1948 //      flowtrack = new AliFlowTrack();
1949 //      flowtrack->SetPhi(fTrackPhi);
1950 //      flowtrack->SetEta(fTrackEta);
1951 //      flowtrack->SetWeight(fTrackWeight);
1952 //      flowtrack->SetPt(fMCparticle->Pt());
1953 //      break;
1954 //    case kTrackWithPtFromFirstMother:
1955 //      if (!fMCparticle) return NULL;
1956 //      flowtrack = new AliFlowTrack();
1957 //      flowtrack->SetPhi(fTrackPhi);
1958 //      flowtrack->SetEta(fTrackEta);
1959 //      flowtrack->SetWeight(fTrackWeight);
1960 //      tmpTParticle = fMCparticle->Particle();
1961 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1962 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
1963 //      break;
1964 //    default:
1965 //      flowtrack = new AliFlowTrack();
1966 //      flowtrack->SetPhi(fTrackPhi);
1967 //      flowtrack->SetEta(fTrackEta);
1968 //      flowtrack->SetWeight(fTrackWeight);
1969 //      break;
1970 //  }
1971 //  flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1972 //  return flowtrack;
1973 //}
1974 //
1975 ////-----------------------------------------------------------------------
1976 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
1977 //{
1978 //  //make flow track from AliVParticle (ESD,AOD,MC)
1979 //  if (!fTrack) return NULL;
1980 //  AliFlowTrack* flowtrack=NULL;
1981 //  TParticle *tmpTParticle=NULL;
1982 //  AliMCParticle* tmpAliMCParticle=NULL;
1983 //  AliExternalTrackParam* externalParams=NULL;
1984 //  AliESDtrack* esdtrack=NULL;
1985 //  switch(fParamMix)
1986 //  {
1987 //    case kPure:
1988 //      flowtrack = new AliFlowTrack(fTrack);
1989 //      break;
1990 //    case kTrackWithMCkine:
1991 //      flowtrack = new AliFlowTrack(fMCparticle);
1992 //      break;
1993 //    case kTrackWithMCPID:
1994 //      flowtrack = new AliFlowTrack(fTrack);
1995 //      //flowtrack->setPID(...) from mc, when implemented
1996 //      break;
1997 //    case kTrackWithMCpt:
1998 //      if (!fMCparticle) return NULL;
1999 //      flowtrack = new AliFlowTrack(fTrack);
2000 //      flowtrack->SetPt(fMCparticle->Pt());
2001 //      break;
2002 //    case kTrackWithPtFromFirstMother:
2003 //      if (!fMCparticle) return NULL;
2004 //      flowtrack = new AliFlowTrack(fTrack);
2005 //      tmpTParticle = fMCparticle->Particle();
2006 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2007 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2008 //      break;
2009 //    case kTrackWithTPCInnerParams:
2010 //      esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2011 //      if (!esdtrack) return NULL;
2012 //      externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2013 //      if (!externalParams) return NULL;
2014 //      flowtrack = new AliFlowTrack(externalParams);
2015 //      break;
2016 //    default:
2017 //      flowtrack = new AliFlowTrack(fTrack);
2018 //      break;
2019 //  }
2020 //  if (fParamType==kMC)
2021 //  {
2022 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
2023 //    flowtrack->SetID(fTrackLabel);
2024 //  }
2025 //  else if (dynamic_cast<AliESDtrack*>(fTrack))
2026 //  {
2027 //    flowtrack->SetSource(AliFlowTrack::kFromESD);
2028 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2029 //  }
2030 //  else if (dynamic_cast<AliAODTrack*>(fTrack))
2031 //  {
2032 //    flowtrack->SetSource(AliFlowTrack::kFromAOD);
2033 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2034 //  }
2035 //  else if (dynamic_cast<AliMCParticle*>(fTrack))
2036 //  {
2037 //    flowtrack->SetSource(AliFlowTrack::kFromMC);
2038 //    flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2039 //  }
2040 //  return flowtrack;
2041 //}
2042 //
2043 ////-----------------------------------------------------------------------
2044 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
2045 //{
2046 //  //make a flow track from PMD track
2047 //  AliFlowTrack* flowtrack=NULL;
2048 //  TParticle *tmpTParticle=NULL;
2049 //  AliMCParticle* tmpAliMCParticle=NULL;
2050 //  switch (fParamMix)
2051 //  {
2052 //    case kPure:
2053 //      flowtrack = new AliFlowTrack();
2054 //      flowtrack->SetPhi(fTrackPhi);
2055 //      flowtrack->SetEta(fTrackEta);
2056 //      flowtrack->SetWeight(fTrackWeight);
2057 //      break;
2058 //    case kTrackWithMCkine:
2059 //      if (!fMCparticle) return NULL;
2060 //      flowtrack = new AliFlowTrack();
2061 //      flowtrack->SetPhi( fMCparticle->Phi() );
2062 //      flowtrack->SetEta( fMCparticle->Eta() );
2063 //      flowtrack->SetWeight(fTrackWeight);
2064 //      flowtrack->SetPt( fMCparticle->Pt() );
2065 //      break;
2066 //    case kTrackWithMCpt:
2067 //      if (!fMCparticle) return NULL;
2068 //      flowtrack = new AliFlowTrack();
2069 //      flowtrack->SetPhi(fTrackPhi);
2070 //      flowtrack->SetEta(fTrackEta);
2071 //      flowtrack->SetWeight(fTrackWeight);
2072 //      flowtrack->SetPt(fMCparticle->Pt());
2073 //      break;
2074 //    case kTrackWithPtFromFirstMother:
2075 //      if (!fMCparticle) return NULL;
2076 //      flowtrack = new AliFlowTrack();
2077 //      flowtrack->SetPhi(fTrackPhi);
2078 //      flowtrack->SetEta(fTrackEta);
2079 //      flowtrack->SetWeight(fTrackWeight);
2080 //      tmpTParticle = fMCparticle->Particle();
2081 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2082 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2083 //      break;
2084 //    default:
2085 //      flowtrack = new AliFlowTrack();
2086 //      flowtrack->SetPhi(fTrackPhi);
2087 //      flowtrack->SetEta(fTrackEta);
2088 //      flowtrack->SetWeight(fTrackWeight);
2089 //      break;
2090 //  }
2091 //
2092 //  flowtrack->SetSource(AliFlowTrack::kFromPMD);
2093 //  return flowtrack;
2094 //}
2095 //
2096 ////-----------------------------------------------------------------------
2097 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVZERO() const
2098 //{
2099 //  //make a flow track from VZERO
2100 //  AliFlowTrack* flowtrack=NULL;
2101 //  TParticle *tmpTParticle=NULL;
2102 //  AliMCParticle* tmpAliMCParticle=NULL;
2103 //  switch (fParamMix)
2104 //  {
2105 //    case kPure:
2106 //      flowtrack = new AliFlowTrack();
2107 //      flowtrack->SetPhi(fTrackPhi);
2108 //      flowtrack->SetEta(fTrackEta);
2109 //      flowtrack->SetWeight(fTrackWeight);
2110 //      break;
2111 //    case kTrackWithMCkine:
2112 //      if (!fMCparticle) return NULL;
2113 //      flowtrack = new AliFlowTrack();
2114 //      flowtrack->SetPhi( fMCparticle->Phi() );
2115 //      flowtrack->SetEta( fMCparticle->Eta() );
2116 //      flowtrack->SetWeight(fTrackWeight);
2117 //      flowtrack->SetPt( fMCparticle->Pt() );
2118 //      break;
2119 //    case kTrackWithMCpt:
2120 //      if (!fMCparticle) return NULL;
2121 //      flowtrack = new AliFlowTrack();
2122 //      flowtrack->SetPhi(fTrackPhi);
2123 //      flowtrack->SetEta(fTrackEta);
2124 //      flowtrack->SetWeight(fTrackWeight);
2125 //      flowtrack->SetPt(fMCparticle->Pt());
2126 //      break;
2127 //    case kTrackWithPtFromFirstMother:
2128 //      if (!fMCparticle) return NULL;
2129 //      flowtrack = new AliFlowTrack();
2130 //      flowtrack->SetPhi(fTrackPhi);
2131 //      flowtrack->SetEta(fTrackEta);
2132 //      flowtrack->SetWeight(fTrackWeight);
2133 //      tmpTParticle = fMCparticle->Particle();
2134 //      tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2135 //      flowtrack->SetPt(tmpAliMCParticle->Pt());
2136 //      break;
2137 //    default:
2138 //      flowtrack = new AliFlowTrack();
2139 //      flowtrack->SetPhi(fTrackPhi);
2140 //      flowtrack->SetEta(fTrackEta);
2141 //      flowtrack->SetWeight(fTrackWeight);
2142 //      break;
2143 //  }
2144 //
2145 //  flowtrack->SetSource(AliFlowTrack::kFromVZERO);
2146 //  return flowtrack;
2147 //}
2148 //
2149 ////-----------------------------------------------------------------------
2150 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
2151 //{
2152 //  //get a flow track constructed from whatever we applied cuts on
2153 //  //caller is resposible for deletion
2154 //  //if construction fails return NULL
2155 //  //TODO: for tracklets, PMD and VZERO we probably need just one method,
2156 //  //something like MakeFlowTrackGeneric(), wait with this until
2157 //  //requirements quirks are known.
2158 //  switch (fParamType)
2159 //  {
2160 //    case kSPDtracklet:
2161 //      return MakeFlowTrackSPDtracklet();
2162 //    case kPMD:
2163 //      return MakeFlowTrackPMDtrack();
2164 //    case kVZERO:
2165 //      return MakeFlowTrackVZERO();
2166 //    default:
2167 //      return MakeFlowTrackVParticle();
2168 //  }
2169 //}
2170
2171 //-----------------------------------------------------------------------
2172 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
2173 {
2174   //check if current particle is a physical primary
2175   if (!fMCevent) return kFALSE;
2176   if (fTrackLabel<0) return kFALSE;
2177   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
2178 }
2179
2180 //-----------------------------------------------------------------------
2181 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
2182 {
2183   //check if current particle is a physical primary
2184   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
2185   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
2186   if (!track) return kFALSE;
2187   TParticle* particle = track->Particle();
2188   Bool_t transported = particle->TestBit(kTransportBit);
2189   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
2190         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
2191   return (physprim && (transported || !requiretransported));
2192 }
2193
2194 //-----------------------------------------------------------------------
2195 void AliFlowTrackCuts::DefineHistograms()
2196 {
2197   //define qa histograms
2198   if (fQA) return;
2199
2200   const Int_t kNbinsP=200;
2201   Double_t binsP[kNbinsP+1];
2202   binsP[0]=0.0;
2203   for(int i=1; i<kNbinsP+1; i++)
2204   {
2205     //if(binsP[i-1]+0.05<1.01)
2206     //  binsP[i]=binsP[i-1]+0.05;
2207     //else
2208     binsP[i]=binsP[i-1]+0.05;
2209   }
2210
2211   const Int_t nBinsDCA=1000;
2212   Double_t binsDCA[nBinsDCA+1];
2213   for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
2214   //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
2215
2216   Bool_t adddirstatus = TH1::AddDirectoryStatus();
2217   TH1::AddDirectory(kFALSE);
2218   fQA=new TList(); fQA->SetOwner();
2219   fQA->SetName(Form("%s QA",GetName()));
2220   TList* before = new TList(); before->SetOwner();
2221   before->SetName("before");
2222   TList* after = new TList(); after->SetOwner();
2223   after->SetName("after");
2224   fQA->Add(before);
2225   fQA->Add(after);
2226   before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2227   after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2228   before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2229   after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2230   before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2231   after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2232   //primary
2233   TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2234   TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2235   TAxis* axis = NULL;
2236   axis = hb->GetYaxis();
2237   axis->SetBinLabel(1,"secondary");
2238   axis->SetBinLabel(2,"primary");
2239   axis = ha->GetYaxis();
2240   axis->SetBinLabel(1,"secondary");
2241   axis->SetBinLabel(2,"primary");
2242   before->Add(hb); //3
2243   after->Add(ha); //3
2244   //production process
2245   hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2246                       -0.5, kMaxMCProcess-0.5);
2247   ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2248                       -0.5, kMaxMCProcess-0.5);
2249   axis = hb->GetYaxis();
2250   for (Int_t i=0; i<kMaxMCProcess; i++)
2251   {
2252     axis->SetBinLabel(i+1,TMCProcessName[i]);
2253   }
2254   axis = ha->GetYaxis();
2255   for (Int_t i=0; i<kMaxMCProcess; i++)
2256   {
2257     axis->SetBinLabel(i+1,TMCProcessName[i]);
2258   }
2259   before->Add(hb); //4
2260   after->Add(ha); //4
2261   //DCA
2262   before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2263   after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2264   before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2265   after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2266   //first mother
2267   hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2268   ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2269   hb->GetYaxis()->SetBinLabel(1, "#gamma");
2270   ha->GetYaxis()->SetBinLabel(1, "#gamma");
2271   hb->GetYaxis()->SetBinLabel(2, "e^{+}");
2272   ha->GetYaxis()->SetBinLabel(2, "e^{+}");
2273   hb->GetYaxis()->SetBinLabel(3, "e^{-}");
2274   ha->GetYaxis()->SetBinLabel(3, "e^{-}");
2275   hb->GetYaxis()->SetBinLabel(4, "#nu");
2276   ha->GetYaxis()->SetBinLabel(4, "#nu");
2277   hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2278   ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2279   hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2280   ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2281   hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2282   ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2283   hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2284   ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2285   hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2286   ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2287   hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2288   ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2289   hb->GetYaxis()->SetBinLabel(11, "K^{+}");
2290   ha->GetYaxis()->SetBinLabel(11, "K^{+}");
2291   hb->GetYaxis()->SetBinLabel(12, "K^{-}");
2292   ha->GetYaxis()->SetBinLabel(12, "K^{-}");
2293   hb->GetYaxis()->SetBinLabel( 13, "n");
2294   ha->GetYaxis()->SetBinLabel( 13, "n");
2295   hb->GetYaxis()->SetBinLabel( 14, "p");
2296   ha->GetYaxis()->SetBinLabel( 14, "p");
2297   hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
2298   ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
2299   hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2300   ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2301   hb->GetYaxis()->SetBinLabel(17, "#eta");
2302   ha->GetYaxis()->SetBinLabel(17, "#eta");
2303   hb->GetYaxis()->SetBinLabel(18, "#Lambda");
2304   ha->GetYaxis()->SetBinLabel(18, "#Lambda");
2305   hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2306   ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2307   hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2308   ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2309   hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2310   ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2311   hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2312   ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2313   hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2314   ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2315   hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2316   ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2317   hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
2318   ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
2319   hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2320   ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2321   hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2322   ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2323   hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2324   ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2325   hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2326   ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2327   hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2328   ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2329   hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2330   ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2331   hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2332   ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2333   hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2334   ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2335   hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2336   ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2337   hb->GetYaxis()->SetBinLabel(35, "D^{+}");
2338   ha->GetYaxis()->SetBinLabel(35, "D^{+}");
2339   hb->GetYaxis()->SetBinLabel(36, "D^{-}");
2340   ha->GetYaxis()->SetBinLabel(36, "D^{-}");
2341   hb->GetYaxis()->SetBinLabel(37, "D^{0}");
2342   ha->GetYaxis()->SetBinLabel(37, "D^{0}");
2343   hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2344   ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2345   hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2346   ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2347   hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2348   ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2349   hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2350   ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2351   hb->GetYaxis()->SetBinLabel(42, "W^{+}");
2352   ha->GetYaxis()->SetBinLabel(42, "W^{+}");
2353   hb->GetYaxis()->SetBinLabel(43, "W^{-}");
2354   ha->GetYaxis()->SetBinLabel(43, "W^{-}");
2355   hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
2356   ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
2357   before->Add(hb);//7
2358   after->Add(ha);//7
2359
2360   before->Add(new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2361   after->Add( new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2362
2363   before->Add(new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
2364   after->Add( new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//9
2365
2366   before->Add(new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
2367   after->Add( new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//10
2368
2369   before->Add(new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
2370   after->Add( new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT",     kNbinsP,binsP,1000,-2e4, 2e4));//11
2371
2372   before->Add(new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
2373   after->Add( new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT",   kNbinsP,binsP,1000,-2e4, 2e4));//12
2374
2375   //kink stuff
2376   before->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2377   after->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2378
2379   before->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2380   after->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2381
2382   before->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2383   after->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2384
2385   before->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2386   after->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2387
2388   before->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2389   after->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2390   TH1::AddDirectory(adddirstatus);
2391 }
2392
2393 //-----------------------------------------------------------------------
2394 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
2395 {
2396   //get the number of tracks in the input event according source
2397   //selection (ESD tracks, tracklets, MC particles etc.)
2398   AliESDEvent* esd=NULL;
2399   AliAODEvent* aod=NULL;  // XZhang 20120615
2400   switch (fParamType)
2401   {
2402     case kSPDtracklet:
2403       if (!fEvent) return 0;                                           // XZhang 20120615
2404       esd = dynamic_cast<AliESDEvent*>(fEvent);
2405       aod = dynamic_cast<AliAODEvent*>(fEvent);                        // XZhang 20120615
2406 //    if (!esd) return 0;                                              // XZhang 20120615
2407 //    return esd->GetMultiplicity()->GetNumberOfTracklets();           // XZhang 20120615
2408       if (esd) return esd->GetMultiplicity()->GetNumberOfTracklets();  // XZhang 20120615
2409       if (aod) return aod->GetTracklets()->GetNumberOfTracklets();     // XZhang 20120615
2410     case kMC:
2411       if (!fMCevent) return 0;
2412       return fMCevent->GetNumberOfTracks();
2413     case kPMD:
2414       esd = dynamic_cast<AliESDEvent*>(fEvent);
2415       if (!esd) return 0;
2416       return esd->GetNumberOfPmdTracks();
2417     case kVZERO:
2418       return fgkNumberOfVZEROtracks;
2419     case kMUON:                                      // XZhang 20120604
2420       if (!fEvent) return 0;                         // XZhang 20120604
2421       esd = dynamic_cast<AliESDEvent*>(fEvent);      // XZhang 20120604
2422       if (esd) return esd->GetNumberOfMuonTracks();  // XZhang 20120604
2423       return fEvent->GetNumberOfTracks();  // if AOD // XZhang 20120604
2424     case kKink:
2425       esd = dynamic_cast<AliESDEvent*>(fEvent);
2426       if (!esd) return 0;
2427       return esd->GetNumberOfKinks();
2428     case kV0:
2429       esd = dynamic_cast<AliESDEvent*>(fEvent);
2430       if (!esd) return 0;
2431       return esd->GetNumberOfV0s();
2432     default:
2433       if (!fEvent) return 0;
2434       return fEvent->GetNumberOfTracks();
2435   }
2436   return 0;
2437 }
2438
2439 //-----------------------------------------------------------------------
2440 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
2441 {
2442   //get the input object according the data source selection:
2443   //(esd tracks, traclets, mc particles,etc...)
2444   AliESDEvent* esd=NULL;
2445   AliAODEvent* aod=NULL;  // XZhang 20120615
2446   switch (fParamType)
2447   {
2448     case kSPDtracklet:
2449       if (!fEvent) return NULL;                                              // XZhang 20120615
2450       esd = dynamic_cast<AliESDEvent*>(fEvent);
2451       aod = dynamic_cast<AliAODEvent*>(fEvent);                              // XZhang 20120615
2452 //    if (!esd) return NULL;                                                 // XZhang 20120615
2453 //    return const_cast<AliMultiplicity*>(esd->GetMultiplicity());           // XZhang 20120615
2454       if (esd) return const_cast<AliMultiplicity*>(esd->GetMultiplicity());  // XZhang 20120615
2455       if (aod) return const_cast<AliAODTracklets*>(aod->GetTracklets());     // XZhang 20120615
2456     case kMC:
2457       if (!fMCevent) return NULL;
2458       return fMCevent->GetTrack(i);
2459     case kPMD:
2460       esd = dynamic_cast<AliESDEvent*>(fEvent);
2461       if (!esd) return NULL;
2462       return esd->GetPmdTrack(i);
2463     case kVZERO:
2464       esd = dynamic_cast<AliESDEvent*>(fEvent);
2465       if (!esd) //contributed by G.Ortona
2466       {
2467         aod = dynamic_cast<AliAODEvent*>(fEvent);
2468         if(!aod)return NULL;
2469         return aod->GetVZEROData();
2470       }
2471       return esd->GetVZEROData();
2472     case kMUON:                                  // XZhang 20120604
2473       if (!fEvent) return NULL;                  // XZhang 20120604
2474       esd = dynamic_cast<AliESDEvent*>(fEvent);  // XZhang 20120604
2475       if (esd) return esd->GetMuonTrack(i);      // XZhang 20120604
2476       return fEvent->GetTrack(i);  // if AOD     // XZhang 20120604
2477     case kKink:
2478       esd = dynamic_cast<AliESDEvent*>(fEvent);
2479       if (!esd) return NULL;
2480       return esd->GetKink(i);
2481     case kV0:
2482       esd = dynamic_cast<AliESDEvent*>(fEvent);
2483       if (!esd) return NULL;
2484       return esd->GetV0(i);
2485     default:
2486       if (!fEvent) return NULL;
2487       return fEvent->GetTrack(i);
2488   }
2489 }
2490
2491 //-----------------------------------------------------------------------
2492 void AliFlowTrackCuts::Clear(Option_t*)
2493 {
2494   //clean up
2495   fMCevent=NULL;
2496   fEvent=NULL;
2497   ClearTrack();
2498 }
2499
2500 //-----------------------------------------------------------------------
2501 void AliFlowTrackCuts::ClearTrack(Option_t*)
2502 {
2503   //clean up last track
2504   fKink=NULL;
2505   fV0=NULL;
2506   fTrack=NULL;
2507   fMCparticle=NULL;
2508   fTrackLabel=-997;
2509   fTrackWeight=1.0;
2510   fTrackEta=0.0;
2511   fTrackPhi=0.0;
2512   fTrackPt=0.0;
2513   fPOItype=1;
2514   fTrackMass=0.;
2515 }
2516 //-----------------------------------------------------------------------
2517 Bool_t AliFlowTrackCuts::PassesAODpidCut(const AliAODTrack* track )
2518 {
2519      if(!track->GetAODEvent()->GetTOFHeader()){
2520           AliAODPid *pidObj = track->GetDetPid();
2521           if (!pidObj) fESDpid.GetTOFResponse().SetTimeResolution(84.);
2522           else{
2523             Double_t sigmaTOFPidInAOD[10];
2524             pidObj->GetTOFpidResolution(sigmaTOFPidInAOD);
2525             if(sigmaTOFPidInAOD[0] > 84.){
2526               fESDpid.GetTOFResponse().SetTimeResolution(sigmaTOFPidInAOD[0]); // use the electron TOF PID sigma as time resolution (including the T0 used)
2527           }
2528         }
2529      }
2530
2531  //check if passes the selected pid cut for ESDs
2532   Bool_t pass = kTRUE;
2533   switch (fPIDsource)
2534   {
2535    case kTOFbeta:
2536       if (!PassesTOFbetaCut(track)) pass=kFALSE;
2537       break;
2538   case kTOFbayesian:
2539       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2540       break;
2541   case kTPCbayesian:
2542       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2543       break;
2544   default:
2545     return kTRUE;
2546     break;
2547  }
2548   return pass;
2549
2550 }
2551 //-----------------------------------------------------------------------
2552 Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
2553 {
2554   //check if passes the selected pid cut for ESDs
2555   Bool_t pass = kTRUE;
2556   switch (fPIDsource)
2557   {
2558     case kTPCpid:
2559       if (!PassesTPCpidCut(track)) pass=kFALSE;
2560       break;
2561     case kTPCdedx:
2562       if (!PassesTPCdedxCut(track)) pass=kFALSE;
2563       break;
2564     case kTOFpid:
2565       if (!PassesTOFpidCut(track)) pass=kFALSE;
2566       break;
2567     case kTOFbeta:
2568       if (!PassesTOFbetaCut(track)) pass=kFALSE;
2569       break;
2570     case kTOFbetaSimple:
2571       if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
2572       break;
2573     case kTPCbayesian:
2574       if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2575       break;
2576       // part added by F. Noferini
2577     case kTOFbayesian:
2578       if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2579       break;
2580       // end part added by F. Noferini
2581
2582       //part added by Natasha
2583     case kTPCNuclei:
2584       if (!PassesNucleiSelection(track)) pass=kFALSE;
2585       break;
2586       //end part added by Natasha
2587
2588     case kTPCTOFNsigma:
2589       if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
2590       break;
2591     default:
2592       printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
2593       pass=kFALSE;
2594       break;
2595   }
2596   return pass;
2597 }
2598
2599 //-----------------------------------------------------------------------
2600 Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
2601 {
2602   //check if passes PID cut using timing in TOF
2603   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2604                      (track->GetTOFsignal() > 12000) &&
2605                      (track->GetTOFsignal() < 100000) &&
2606                      (track->GetIntegratedLength() > 365);
2607
2608   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2609
2610   Bool_t statusMatchingHard = TPCTOFagree(track);
2611   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2612        return kFALSE;
2613
2614   if (!goodtrack) return kFALSE;
2615
2616   const Float_t c = 2.99792457999999984e-02;
2617   Float_t p = track->GetP();
2618   Float_t l = track->GetIntegratedLength();
2619   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2620   Float_t timeTOF = track->GetTOFsignal()- trackT0;
2621   Float_t beta = l/timeTOF/c;
2622   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2623   track->GetIntegratedTimes(integratedTimes);
2624   Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
2625   Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
2626   for (Int_t i=0;i<5;i++)
2627   {
2628     betaHypothesis[i] = l/integratedTimes[i]/c;
2629     s[i] = beta-betaHypothesis[i];
2630   }
2631
2632   switch (fParticleID)
2633   {
2634     case AliPID::kPion:
2635       return ( (s[2]<0.015) && (s[2]>-0.015) &&
2636                (s[3]>0.025) &&
2637                (s[4]>0.03) );
2638     case AliPID::kKaon:
2639       return ( (s[3]<0.015) && (s[3]>-0.015) &&
2640                (s[2]<-0.03) &&
2641                (s[4]>0.03) );
2642     case AliPID::kProton:
2643       return ( (s[4]<0.015) && (s[4]>-0.015) &&
2644                (s[3]<-0.025) &&
2645                (s[2]<-0.025) );
2646     default:
2647       return kFALSE;
2648   }
2649   return kFALSE;
2650 }
2651
2652 //-----------------------------------------------------------------------
2653 Float_t AliFlowTrackCuts::GetBeta(const AliVTrack* track)
2654 {
2655   //get beta
2656   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2657   track->GetIntegratedTimes(integratedTimes);
2658
2659   const Float_t c = 2.99792457999999984e-02;
2660   Float_t p = track->P();
2661   Float_t l = integratedTimes[0]*c;
2662   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2663   Float_t timeTOF = track->GetTOFsignal()- trackT0;
2664   return l/timeTOF/c;
2665 }
2666 //-----------------------------------------------------------------------
2667 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliAODTrack* track )
2668 {
2669   //check if track passes pid selection with an asymmetric TOF beta cut
2670   if (!fTOFpidCuts)
2671   {
2672     //printf("no TOFpidCuts\n");
2673     return kFALSE;
2674   }
2675
2676   //check if passes PID cut using timing in TOF
2677   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2678                      (track->GetTOFsignal() > 12000) &&
2679                      (track->GetTOFsignal() < 100000);
2680
2681   if (!goodtrack) return kFALSE;
2682
2683   const Float_t c = 2.99792457999999984e-02;
2684   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2685   track->GetIntegratedTimes(integratedTimes);
2686   Float_t l = integratedTimes[0]*c;
2687
2688   goodtrack = goodtrack && (l > 365);
2689
2690   if (!goodtrack) return kFALSE;
2691
2692   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2693
2694   Bool_t statusMatchingHard = TPCTOFagree(track);
2695   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2696        return kFALSE;
2697
2698
2699   Float_t beta = GetBeta(track);
2700
2701   //construct the pid index because it's not AliPID::EParticleType
2702   Int_t pid = 0;
2703   switch (fParticleID)
2704   {
2705     case AliPID::kPion:
2706       pid=2;
2707       break;
2708     case AliPID::kKaon:
2709       pid=3;
2710       break;
2711     case AliPID::kProton:
2712       pid=4;
2713       break;
2714     default:
2715       return kFALSE;
2716   }
2717
2718   //signal to cut on
2719   Float_t p = track->P();
2720   Float_t betahypothesis = l/integratedTimes[pid]/c;
2721   Float_t betadiff = beta-betahypothesis;
2722
2723   Float_t* arr = fTOFpidCuts->GetMatrixArray();
2724   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2725   if (col<0) return kFALSE;
2726   Float_t min = (*fTOFpidCuts)(1,col);
2727   Float_t max = (*fTOFpidCuts)(2,col);
2728
2729   Bool_t pass = (betadiff>min && betadiff<max);
2730
2731   return pass;
2732 }
2733 //-----------------------------------------------------------------------
2734 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
2735 {
2736   //check if track passes pid selection with an asymmetric TOF beta cut
2737   if (!fTOFpidCuts)
2738   {
2739     //printf("no TOFpidCuts\n");
2740     return kFALSE;
2741   }
2742
2743   //check if passes PID cut using timing in TOF
2744   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2745                      (track->GetTOFsignal() > 12000) &&
2746                      (track->GetTOFsignal() < 100000) &&
2747                      (track->GetIntegratedLength() > 365);
2748
2749   if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2750
2751   if (!goodtrack) return kFALSE;
2752
2753   Bool_t statusMatchingHard = TPCTOFagree(track);
2754   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2755        return kFALSE;
2756
2757   Float_t beta = GetBeta(track);
2758   Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
2759   track->GetIntegratedTimes(integratedTimes);
2760
2761   //construct the pid index because it's not AliPID::EParticleType
2762   Int_t pid = 0;
2763   switch (fParticleID)
2764   {
2765     case AliPID::kPion:
2766       pid=2;
2767       break;
2768     case AliPID::kKaon:
2769       pid=3;
2770       break;
2771     case AliPID::kProton:
2772       pid=4;
2773       break;
2774     default:
2775       return kFALSE;
2776   }
2777
2778   //signal to cut on
2779   const Float_t c = 2.99792457999999984e-02;
2780   Float_t l = track->GetIntegratedLength();
2781   Float_t p = track->GetP();
2782   Float_t betahypothesis = l/integratedTimes[pid]/c;
2783   Float_t betadiff = beta-betahypothesis;
2784
2785   Float_t* arr = fTOFpidCuts->GetMatrixArray();
2786   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2787   if (col<0) return kFALSE;
2788   Float_t min = (*fTOFpidCuts)(1,col);
2789   Float_t max = (*fTOFpidCuts)(2,col);
2790
2791   Bool_t pass = (betadiff>min && betadiff<max);
2792
2793   return pass;
2794 }
2795
2796 //-----------------------------------------------------------------------
2797 Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
2798 {
2799   //check if passes PID cut using default TOF pid
2800   Double_t pidTOF[AliPID::kSPECIES];
2801   track->GetTOFpid(pidTOF);
2802   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
2803   return kFALSE;
2804 }
2805
2806 //-----------------------------------------------------------------------
2807 Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
2808 {
2809   //check if passes PID cut using default TPC pid
2810   Double_t pidTPC[AliPID::kSPECIES];
2811   track->GetTPCpid(pidTPC);
2812   Double_t probablity = 0.;
2813   switch (fParticleID)
2814   {
2815     case AliPID::kPion:
2816       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
2817       break;
2818     default:
2819       probablity = pidTPC[fParticleID];
2820   }
2821   if (probablity >= fParticleProbability) return kTRUE;
2822   return kFALSE;
2823 }
2824
2825 //-----------------------------------------------------------------------
2826 Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
2827 {
2828   //get TPC dedx
2829   return track->GetTPCsignal();
2830 }
2831
2832 //-----------------------------------------------------------------------
2833 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
2834 {
2835   //check if passes PID cut using dedx signal in the TPC
2836   if (!fTPCpidCuts)
2837   {
2838     //printf("no TPCpidCuts\n");
2839     return kFALSE;
2840   }
2841
2842   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
2843   if (!tpcparam) return kFALSE;
2844   Double_t p = tpcparam->GetP();
2845   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
2846   Float_t sigTPC = track->GetTPCsignal();
2847   Float_t s = (sigTPC-sigExp)/sigExp;
2848
2849   Float_t* arr = fTPCpidCuts->GetMatrixArray();
2850   Int_t arrSize = fTPCpidCuts->GetNcols();
2851   Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
2852   if (col<0) return kFALSE;
2853   Float_t min = (*fTPCpidCuts)(1,col);
2854   Float_t max = (*fTPCpidCuts)(2,col);
2855
2856   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
2857   return (s>min && s<max);
2858 }
2859
2860 //-----------------------------------------------------------------------
2861 void AliFlowTrackCuts::InitPIDcuts()
2862 {
2863   //init matrices with PID cuts
2864   TMatrixF* t = NULL;
2865   if (!fTPCpidCuts)
2866   {
2867     if (fParticleID==AliPID::kPion)
2868     {
2869       t = new TMatrixF(3,15);
2870       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.0;
2871       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.1;
2872       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.2;
2873       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.2;
2874       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
2875       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
2876       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
2877       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
2878       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
2879       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.4;  (*t)(2,9)  =  0.05;
2880       (*t)(0,10)  = 0.70;  (*t)(1,10)  = -0.4;  (*t)(2,10)  =     0;
2881       (*t)(0,11)  = 0.75;  (*t)(1,11)  = -0.4;  (*t)(2,11)  =     0;
2882       (*t)(0,12)  = 0.80;  (*t)(1,12)  = -0.4;  (*t)(2,12)  = -0.05;
2883       (*t)(0,13)  = 0.85;  (*t)(1,13)  = -0.4;  (*t)(2,13)  = -0.1;
2884       (*t)(0,14)  = 0.90;  (*t)(1,14)  = 0;     (*t)(2,14)  =     0;
2885     }
2886     else
2887     if (fParticleID==AliPID::kKaon)
2888     {
2889       t = new TMatrixF(3,12);
2890       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.2;
2891       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  = 0.2;
2892       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.2;  (*t)(2,2)  = 0.2;
2893       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.2;  (*t)(2,3)  = 0.2;
2894       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.2;
2895       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.2;
2896       (*t)(0,6)  = 0.50;  (*t)(1,6)  =-0.05;  (*t)(2,6)  = 0.2;
2897       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.1;  (*t)(2,7)  = 0.1;
2898       (*t)(0,8)  = 0.60;  (*t)(1,8)  =-0.05;  (*t)(2,8)  = 0.1;
2899       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  = 0.15;
2900       (*t)(0,10)  = 0.70;  (*t)(1,10)  = 0.05;  (*t)(2,10)  = 0.2;
2901       (*t)(0,11)  = 0.75;  (*t)(1,11)  =    0;  (*t)(2,11)  = 0;
2902     }
2903     else
2904     if (fParticleID==AliPID::kProton)
2905     {
2906       t = new TMatrixF(3,9);
2907       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.1;  (*t)(2,0)  =  0.1;
2908       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.2;  (*t)(2,1)  =  0.2;
2909       (*t)(0,2)  = 0.80;  (*t)(1,2)  = -0.1;  (*t)(2,2)  =  0.2;
2910       (*t)(0,3)  = 0.85;  (*t)(1,3)  =-0.05;  (*t)(2,3)  =  0.2;
2911       (*t)(0,4)  = 0.90;  (*t)(1,4)  =-0.05;  (*t)(2,4)  = 0.25;
2912       (*t)(0,5)  = 0.95;  (*t)(1,5)  =-0.05;  (*t)(2,5)  = 0.25;
2913       (*t)(0,6)  = 1.00;  (*t)(1,6)  = -0.1;  (*t)(2,6)  = 0.25;
2914       (*t)(0,7)  = 1.10;  (*t)(1,7)  =-0.05;  (*t)(2,7)  =  0.3;
2915       (*t)(0,8) = 1.20;   (*t)(1,8)  =    0;  (*t)(2,8) =    0;
2916     }
2917     delete fTPCpidCuts;
2918     fTPCpidCuts=t;
2919   }
2920   t = NULL;
2921   if (!fTOFpidCuts)
2922   {
2923     if (fParticleID==AliPID::kPion)
2924     {
2925       //TOF pions, 0.9 purity
2926       t = new TMatrixF(3,61);
2927       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2928       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2929       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2930       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2931       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.030;  (*t)(2,4)  =   0.030;
2932       (*t)(0,5)  = 0.250;  (*t)(1,5)  = -0.036;  (*t)(2,5)  =   0.032;
2933       (*t)(0,6)  = 0.300;  (*t)(1,6)  = -0.038;  (*t)(2,6)  =   0.032;
2934       (*t)(0,7)  = 0.350;  (*t)(1,7)  = -0.034;  (*t)(2,7)  =   0.032;
2935       (*t)(0,8)  = 0.400;  (*t)(1,8)  = -0.032;  (*t)(2,8)  =   0.020;
2936       (*t)(0,9)  = 0.450;  (*t)(1,9)  = -0.030;  (*t)(2,9)  =   0.020;
2937       (*t)(0,10)  = 0.500;  (*t)(1,10)  = -0.030;  (*t)(2,10)  =   0.020;
2938       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.030;  (*t)(2,11)  =   0.020;
2939       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.030;  (*t)(2,12)  =   0.020;
2940       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.030;  (*t)(2,13)  =   0.020;
2941       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.030;  (*t)(2,14)  =   0.020;
2942       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.030;  (*t)(2,15)  =   0.020;
2943       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.030;  (*t)(2,16)  =   0.020;
2944       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.030;  (*t)(2,17)  =   0.020;
2945       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.030;  (*t)(2,18)  =   0.020;
2946       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.028;  (*t)(2,19)  =   0.028;
2947       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.028;  (*t)(2,20)  =   0.028;
2948       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.028;  (*t)(2,21)  =   0.028;
2949       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.028;
2950       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.028;
2951       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.020;  (*t)(2,24)  =   0.028;
2952       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.018;  (*t)(2,25)  =   0.028;
2953       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.016;  (*t)(2,26)  =   0.028;
2954       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.014;  (*t)(2,27)  =   0.028;
2955       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.012;  (*t)(2,28)  =   0.026;
2956       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.010;  (*t)(2,29)  =   0.026;
2957       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.008;  (*t)(2,30)  =   0.026;
2958       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.008;  (*t)(2,31)  =   0.024;
2959       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.006;  (*t)(2,32)  =   0.024;
2960       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.004;  (*t)(2,33)  =   0.024;
2961       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.004;  (*t)(2,34)  =   0.024;
2962       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.002;  (*t)(2,35)  =   0.024;
2963       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.002;  (*t)(2,36)  =   0.024;
2964       (*t)(0,37)  = 2.700;  (*t)(1,37)  = 0.000;  (*t)(2,37)  =   0.024;
2965       (*t)(0,38)  = 2.800;  (*t)(1,38)  = 0.000;  (*t)(2,38)  =   0.026;
2966       (*t)(0,39)  = 2.900;  (*t)(1,39)  = 0.000;  (*t)(2,39)  =   0.024;
2967       (*t)(0,40)  = 3.000;  (*t)(1,40)  = 0.002;  (*t)(2,40)  =   0.026;
2968       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.002;  (*t)(2,41)  =   0.026;
2969       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.002;  (*t)(2,42)  =   0.026;
2970       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.002;  (*t)(2,43)  =   0.026;
2971       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.002;  (*t)(2,44)  =   0.026;
2972       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.002;  (*t)(2,45)  =   0.026;
2973       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.002;  (*t)(2,46)  =   0.026;
2974       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.002;  (*t)(2,47)  =   0.026;
2975       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.002;  (*t)(2,48)  =   0.026;
2976       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.004;  (*t)(2,49)  =   0.024;
2977       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.004;  (*t)(2,50)  =   0.026;
2978       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.004;  (*t)(2,51)  =   0.026;
2979       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.004;  (*t)(2,52)  =   0.024;
2980       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.006;  (*t)(2,53)  =   0.024;
2981       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
2982       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
2983       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
2984       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
2985       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
2986       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
2987       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
2988     }
2989     else
2990     if (fParticleID==AliPID::kProton)
2991     {
2992       //TOF protons, 0.9 purity
2993       t = new TMatrixF(3,61);
2994       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
2995       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
2996       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
2997       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
2998       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.07;  (*t)(2,4)  =   0.07;
2999       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.07;  (*t)(2,5)  =   0.07;
3000       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.07;  (*t)(2,6)  =   0.07;
3001       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.07;  (*t)(2,7)  =   0.07;
3002       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.07;  (*t)(2,8)  =   0.07;
3003       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.07;  (*t)(2,9)  =   0.07;
3004       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.07;  (*t)(2,10)  =   0.07;
3005       (*t)(0,11)  = 0.200;  (*t)(1,11)  = -0.07;  (*t)(2,11)  =   0.07;
3006       (*t)(0,12)  = 0.200;  (*t)(1,12)  = -0.07;  (*t)(2,12)  =   0.07;
3007       (*t)(0,13)  = 0.200;  (*t)(1,13)  = -0.07;  (*t)(2,13)  =   0.07;
3008       (*t)(0,14)  = 0.200;  (*t)(1,14)  = -0.07;  (*t)(2,14)  =   0.07;
3009       (*t)(0,15)  = 0.200;  (*t)(1,15)  = -0.07;  (*t)(2,15)  =   0.07;
3010       (*t)(0,16)  = 0.200;  (*t)(1,16)  = -0.07;  (*t)(2,16)  =   0.07;
3011       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.070;  (*t)(2,17)  =   0.070;
3012       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.072;  (*t)(2,18)  =   0.072;
3013       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.072;  (*t)(2,19)  =   0.072;
3014       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.074;  (*t)(2,20)  =   0.074;
3015       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.032;  (*t)(2,21)  =   0.032;
3016       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.026;  (*t)(2,22)  =   0.026;
3017       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.026;  (*t)(2,23)  =   0.026;
3018       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.024;  (*t)(2,24)  =   0.024;
3019       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.024;  (*t)(2,25)  =   0.024;
3020       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.026;  (*t)(2,26)  =   0.026;
3021       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.026;  (*t)(2,27)  =   0.026;
3022       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.026;  (*t)(2,28)  =   0.026;
3023       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.026;  (*t)(2,29)  =   0.026;
3024       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.026;  (*t)(2,30)  =   0.026;
3025       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.026;
3026       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.026;  (*t)(2,32)  =   0.024;
3027       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.028;  (*t)(2,33)  =   0.022;
3028       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.028;  (*t)(2,34)  =   0.020;
3029       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.028;  (*t)(2,35)  =   0.018;
3030       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.028;  (*t)(2,36)  =   0.016;
3031       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.028;  (*t)(2,37)  =   0.016;
3032       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.030;  (*t)(2,38)  =   0.014;
3033       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.030;  (*t)(2,39)  =   0.012;
3034       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.030;  (*t)(2,40)  =   0.012;
3035       (*t)(0,41)  = 3.100;  (*t)(1,41)  = -0.030;  (*t)(2,41)  =   0.010;
3036       (*t)(0,42)  = 3.200;  (*t)(1,42)  = -0.030;  (*t)(2,42)  =   0.010;
3037       (*t)(0,43)  = 3.300;  (*t)(1,43)  = -0.030;  (*t)(2,43)  =   0.010;
3038       (*t)(0,44)  = 3.400;  (*t)(1,44)  = -0.030;  (*t)(2,44)  =   0.008;
3039       (*t)(0,45)  = 3.500;  (*t)(1,45)  = -0.030;  (*t)(2,45)  =   0.008;
3040       (*t)(0,46)  = 3.600;  (*t)(1,46)  = -0.030;  (*t)(2,46)  =   0.008;
3041       (*t)(0,47)  = 3.700;  (*t)(1,47)  = -0.030;  (*t)(2,47)  =   0.006;
3042       (*t)(0,48)  = 3.800;  (*t)(1,48)  = -0.030;  (*t)(2,48)  =   0.006;
3043       (*t)(0,49)  = 3.900;  (*t)(1,49)  = -0.030;  (*t)(2,49)  =   0.006;
3044       (*t)(0,50)  = 4.000;  (*t)(1,50)  = -0.028;  (*t)(2,50)  =   0.004;
3045       (*t)(0,51)  = 4.100;  (*t)(1,51)  = -0.030;  (*t)(2,51)  =   0.004;
3046       (*t)(0,52)  = 4.200;  (*t)(1,52)  = -0.030;  (*t)(2,52)  =   0.004;
3047       (*t)(0,53)  = 4.300;  (*t)(1,53)  = -0.028;  (*t)(2,53)  =   0.002;
3048       (*t)(0,54)  = 4.400;  (*t)(1,54)  = -0.030;  (*t)(2,54)  =   0.002;
3049       (*t)(0,55)  = 4.500;  (*t)(1,55)  = -0.028;  (*t)(2,55)  =   0.002;
3050       (*t)(0,56)  = 4.600;  (*t)(1,56)  = -0.028;  (*t)(2,56)  =   0.002;
3051       (*t)(0,57)  = 4.700;  (*t)(1,57)  = -0.028;  (*t)(2,57)  =   0.000;
3052       (*t)(0,58)  = 4.800;  (*t)(1,58)  = -0.028;  (*t)(2,58)  =   0.002;
3053       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
3054       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
3055     }
3056     else
3057     if (fParticleID==AliPID::kKaon)
3058     {
3059       //TOF kaons, 0.9 purity
3060       t = new TMatrixF(3,61);
3061       (*t)(0,0)  = 0.000;  (*t)(1,0)  = 0.000;  (*t)(2,0)  =   0.000;
3062       (*t)(0,1)  = 0.050;  (*t)(1,1)  = 0.000;  (*t)(2,1)  =   0.000;
3063       (*t)(0,2)  = 0.100;  (*t)(1,2)  = 0.000;  (*t)(2,2)  =   0.000;
3064       (*t)(0,3)  = 0.150;  (*t)(1,3)  = 0.000;  (*t)(2,3)  =   0.000;
3065       (*t)(0,4)  = 0.200;  (*t)(1,4)  = -0.05;  (*t)(2,4)  =   0.05;
3066       (*t)(0,5)  = 0.200;  (*t)(1,5)  = -0.05;  (*t)(2,5)  =   0.05;
3067       (*t)(0,6)  = 0.200;  (*t)(1,6)  = -0.05;  (*t)(2,6)  =   0.05;
3068       (*t)(0,7)  = 0.200;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =   0.05;
3069       (*t)(0,8)  = 0.200;  (*t)(1,8)  = -0.05;  (*t)(2,8)  =   0.05;
3070       (*t)(0,9)  = 0.200;  (*t)(1,9)  = -0.05;  (*t)(2,9)  =   0.05;
3071       (*t)(0,10)  = 0.200;  (*t)(1,10)  = -0.05;  (*t)(2,10)  =   0.05;
3072       (*t)(0,11)  = 0.550;  (*t)(1,11)  = -0.026;  (*t)(2,11)  =   0.026;
3073       (*t)(0,12)  = 0.600;  (*t)(1,12)  = -0.026;  (*t)(2,12)  =   0.026;
3074       (*t)(0,13)  = 0.650;  (*t)(1,13)  = -0.026;  (*t)(2,13)  =   0.026;
3075       (*t)(0,14)  = 0.700;  (*t)(1,14)  = -0.026;  (*t)(2,14)  =   0.026;
3076       (*t)(0,15)  = 0.750;  (*t)(1,15)  = -0.026;  (*t)(2,15)  =   0.026;
3077       (*t)(0,16)  = 0.800;  (*t)(1,16)  = -0.026;  (*t)(2,16)  =   0.026;
3078       (*t)(0,17)  = 0.850;  (*t)(1,17)  = -0.024;  (*t)(2,17)  =   0.024;
3079       (*t)(0,18)  = 0.900;  (*t)(1,18)  = -0.024;  (*t)(2,18)  =   0.024;
3080       (*t)(0,19)  = 0.950;  (*t)(1,19)  = -0.024;  (*t)(2,19)  =   0.024;
3081       (*t)(0,20)  = 1.000;  (*t)(1,20)  = -0.024;  (*t)(2,20)  =   0.024;
3082       (*t)(0,21)  = 1.100;  (*t)(1,21)  = -0.024;  (*t)(2,21)  =   0.024;
3083       (*t)(0,22)  = 1.200;  (*t)(1,22)  = -0.024;  (*t)(2,22)  =   0.022;
3084       (*t)(0,23)  = 1.300;  (*t)(1,23)  = -0.024;  (*t)(2,23)  =   0.020;
3085       (*t)(0,24)  = 1.400;  (*t)(1,24)  = -0.026;  (*t)(2,24)  =   0.016;
3086       (*t)(0,25)  = 1.500;  (*t)(1,25)  = -0.028;  (*t)(2,25)  =   0.014;
3087       (*t)(0,26)  = 1.600;  (*t)(1,26)  = -0.028;  (*t)(2,26)  =   0.012;
3088       (*t)(0,27)  = 1.700;  (*t)(1,27)  = -0.028;  (*t)(2,27)  =   0.010;
3089       (*t)(0,28)  = 1.800;  (*t)(1,28)  = -0.028;  (*t)(2,28)  =   0.010;
3090       (*t)(0,29)  = 1.900;  (*t)(1,29)  = -0.028;  (*t)(2,29)  =   0.008;
3091       (*t)(0,30)  = 2.000;  (*t)(1,30)  = -0.028;  (*t)(2,30)  =   0.006;
3092       (*t)(0,31)  = 2.100;  (*t)(1,31)  = -0.026;  (*t)(2,31)  =   0.006;
3093       (*t)(0,32)  = 2.200;  (*t)(1,32)  = -0.024;  (*t)(2,32)  =   0.004;
3094       (*t)(0,33)  = 2.300;  (*t)(1,33)  = -0.020;  (*t)(2,33)  =   0.002;
3095       (*t)(0,34)  = 2.400;  (*t)(1,34)  = -0.020;  (*t)(2,34)  =   0.002;
3096       (*t)(0,35)  = 2.500;  (*t)(1,35)  = -0.018;  (*t)(2,35)  =   0.000;
3097       (*t)(0,36)  = 2.600;  (*t)(1,36)  = -0.016;  (*t)(2,36)  =   0.000;
3098       (*t)(0,37)  = 2.700;  (*t)(1,37)  = -0.014;  (*t)(2,37)  =   -0.002;
3099       (*t)(0,38)  = 2.800;  (*t)(1,38)  = -0.014;  (*t)(2,38)  =   -0.004;
3100       (*t)(0,39)  = 2.900;  (*t)(1,39)  = -0.012;  (*t)(2,39)  =   -0.004;
3101       (*t)(0,40)  = 3.000;  (*t)(1,40)  = -0.010;  (*t)(2,40)  =   -0.006;
3102       (*t)(0,41)  = 3.100;  (*t)(1,41)  = 0.000;  (*t)(2,41)  =   0.000;
3103       (*t)(0,42)  = 3.200;  (*t)(1,42)  = 0.000;  (*t)(2,42)  =   0.000;
3104       (*t)(0,43)  = 3.300;  (*t)(1,43)  = 0.000;  (*t)(2,43)  =   0.000;
3105       (*t)(0,44)  = 3.400;  (*t)(1,44)  = 0.000;  (*t)(2,44)  =   0.000;
3106       (*t)(0,45)  = 3.500;  (*t)(1,45)  = 0.000;  (*t)(2,45)  =   0.000;
3107       (*t)(0,46)  = 3.600;  (*t)(1,46)  = 0.000;  (*t)(2,46)  =   0.000;
3108       (*t)(0,47)  = 3.700;  (*t)(1,47)  = 0.000;  (*t)(2,47)  =   0.000;
3109       (*t)(0,48)  = 3.800;  (*t)(1,48)  = 0.000;  (*t)(2,48)  =   0.000;
3110       (*t)(0,49)  = 3.900;  (*t)(1,49)  = 0.000;  (*t)(2,49)  =   0.000;
3111       (*t)(0,50)  = 4.000;  (*t)(1,50)  = 0.000;  (*t)(2,50)  =   0.000;
3112       (*t)(0,51)  = 4.100;  (*t)(1,51)  = 0.000;  (*t)(2,51)  =   0.000;
3113       (*t)(0,52)  = 4.200;  (*t)(1,52)  = 0.000;  (*t)(2,52)  =   0.000;
3114       (*t)(0,53)  = 4.300;  (*t)(1,53)  = 0.000;  (*t)(2,53)  =   0.000;
3115       (*t)(0,54)  = 4.400;  (*t)(1,54)  = 0.000;  (*t)(2,54)  =   0.000;
3116       (*t)(0,55)  = 4.500;  (*t)(1,55)  = 0.000;  (*t)(2,55)  =   0.000;
3117       (*t)(0,56)  = 4.600;  (*t)(1,56)  = 0.000;  (*t)(2,56)  =   0.000;
3118       (*t)(0,57)  = 4.700;  (*t)(1,57)  = 0.000;  (*t)(2,57)  =   0.000;
3119       (*t)(0,58)  = 4.800;  (*t)(1,58)  = 0.000;  (*t)(2,58)  =   0.000;
3120       (*t)(0,59)  = 4.900;  (*t)(1,59)  = 0.000;  (*t)(2,59)  =   0.000;
3121       (*t)(0,60)  = 5.900;  (*t)(1,60)  = 0.000;  (*t)(2,60)  =   0.000;
3122     }
3123     delete fTOFpidCuts;
3124     fTOFpidCuts=t;
3125   }
3126 }
3127
3128 //-----------------------------------------------------------------------
3129 //-----------------------------------------------------------------------
3130 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliAODTrack* track)
3131 {
3132   fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3133   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3134
3135   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3136
3137   if(! kTPC) return kFALSE;
3138
3139   fProbBayes = 0.0;
3140
3141 switch (fParticleID)
3142   {
3143     case AliPID::kPion:
3144       fProbBayes = probabilities[2];
3145       break;
3146     case AliPID::kKaon:
3147        fProbBayes = probabilities[3];
3148      break;
3149     case AliPID::kProton:
3150        fProbBayes = probabilities[4];
3151       break;
3152     case AliPID::kElectron:
3153        fProbBayes = probabilities[0];
3154      break;
3155     case AliPID::kMuon:
3156        fProbBayes = probabilities[1];
3157      break;
3158     case AliPID::kDeuteron:
3159        fProbBayes = probabilities[5];
3160      break;
3161     case AliPID::kTriton:
3162        fProbBayes = probabilities[6];
3163      break;
3164     case AliPID::kHe3:
3165        fProbBayes = probabilities[7];
3166      break;
3167    default:
3168       return kFALSE;
3169   }
3170
3171  if(fProbBayes > fParticleProbability){
3172     if(!fCutCharge)
3173       return kTRUE;
3174     else if (fCutCharge && fCharge * track->Charge() > 0)
3175       return kTRUE;
3176   }
3177   return kFALSE;
3178 }
3179 //-----------------------------------------------------------------------
3180 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
3181 {
3182   //cut on TPC bayesian pid
3183   //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3184
3185   //Bool_t statusMatchingHard = TPCTOFagree(track);
3186   //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3187   //     return kFALSE;
3188   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3189   Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3190   //Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
3191
3192   if(! kTPC) return kFALSE;
3193
3194   //  Bool_t statusMatchingHard = 1;
3195   //  Float_t mismProb = 0;
3196   //  if(kTOF){
3197   //    statusMatchingHard = TPCTOFagree(track);
3198   //    mismProb = fBayesianResponse->GetTOFMismProb();
3199   //  }
3200   //  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3201   //       return kFALSE;
3202
3203   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3204
3205   fProbBayes = 0.0;
3206
3207   switch (fParticleID)
3208   {
3209     case AliPID::kPion:
3210       fProbBayes = probabilities[2];
3211       break;
3212     case AliPID::kKaon:
3213       fProbBayes = probabilities[3];
3214      break;
3215     case AliPID::kProton:
3216       fProbBayes = probabilities[4];
3217       break;
3218     case AliPID::kElectron:
3219        fProbBayes = probabilities[0];
3220      break;
3221     case AliPID::kMuon:
3222        fProbBayes = probabilities[1];
3223      break;
3224     case AliPID::kDeuteron:
3225        fProbBayes = probabilities[5];
3226      break;
3227     case AliPID::kTriton:
3228        fProbBayes = probabilities[6];
3229      break;
3230     case AliPID::kHe3:
3231        fProbBayes = probabilities[7];
3232      break;
3233     default:
3234       return kFALSE;
3235   }
3236
3237   if(fProbBayes > fParticleProbability)
3238     {
3239       if(!fCutCharge)
3240         return kTRUE;
3241       else if (fCutCharge && fCharge * track->GetSign() > 0)
3242         return kTRUE;
3243     }
3244   return kFALSE;
3245 }
3246 //-----------------------------------------------------------------------
3247 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliAODTrack* track)
3248 {
3249 //check is track passes bayesian combined TOF+TPC pid cut
3250   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3251                      (track->GetStatus() & AliESDtrack::kTIME) &&
3252                      (track->GetTOFsignal() > 12000) &&
3253                      (track->GetTOFsignal() < 100000);
3254
3255   if (! goodtrack)
3256        return kFALSE;
3257
3258   Bool_t statusMatchingHard = TPCTOFagree(track);
3259 //ciao
3260   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3261        return kFALSE;
3262
3263  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3264   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3265
3266   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3267
3268   fProbBayes = 0.0;
3269
3270  switch (fParticleID)
3271   {
3272     case AliPID::kPion:
3273       fProbBayes = probabilities[2];
3274       break;
3275     case AliPID::kKaon:
3276        fProbBayes = probabilities[3];
3277      break;
3278     case AliPID::kProton:
3279        fProbBayes = probabilities[4];
3280       break;
3281     case AliPID::kElectron:
3282        fProbBayes = probabilities[0];
3283      break;
3284     case AliPID::kMuon:
3285        fProbBayes = probabilities[1];
3286      break;
3287     case AliPID::kDeuteron:
3288        fProbBayes = probabilities[5];
3289      break;
3290     case AliPID::kTriton:
3291        fProbBayes = probabilities[6];
3292      break;
3293     case AliPID::kHe3:
3294        fProbBayes = probabilities[7];
3295      break;
3296    default:
3297       return kFALSE;
3298   }
3299
3300  if(fProbBayes > fParticleProbability && mismProb < 0.5){
3301     if(!fCutCharge)
3302       return kTRUE;
3303     else if (fCutCharge && fCharge * track->Charge() > 0)
3304       return kTRUE;
3305   }
3306   return kFALSE;
3307
3308 }
3309 //-----------------------------------------------------------------------
3310 // part added by F. Noferini (some methods)
3311 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
3312 {
3313   //check is track passes bayesian combined TOF+TPC pid cut
3314   Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3315                      (track->GetStatus() & AliESDtrack::kTIME) &&
3316                      (track->GetTOFsignal() > 12000) &&
3317                      (track->GetTOFsignal() < 100000) &&
3318                      (track->GetIntegratedLength() > 365);
3319
3320   if (! goodtrack)
3321        return kFALSE;
3322
3323   if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3324
3325   Bool_t statusMatchingHard = TPCTOFagree(track);
3326   if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3327        return kFALSE;
3328
3329   fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3330   Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3331
3332   Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3333
3334   fProbBayes = 0.0;
3335
3336   switch (fParticleID)
3337   {
3338     case AliPID::kPion:
3339       fProbBayes = probabilities[2];
3340       break;
3341     case AliPID::kKaon:
3342        fProbBayes = probabilities[3];
3343      break;
3344     case AliPID::kProton:
3345        fProbBayes = probabilities[4];
3346       break;
3347     case AliPID::kElectron:
3348        fProbBayes = probabilities[0];
3349      break;
3350     case AliPID::kMuon:
3351        fProbBayes = probabilities[1];
3352      break;
3353     case AliPID::kDeuteron:
3354        fProbBayes = probabilities[5];
3355      break;
3356     case AliPID::kTriton:
3357        fProbBayes = probabilities[6];
3358      break;
3359     case AliPID::kHe3:
3360        fProbBayes = probabilities[7];
3361      break;
3362    default:
3363       return kFALSE;
3364   }
3365
3366   //  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);
3367   if(fProbBayes > fParticleProbability && mismProb < 0.5){
3368     if(!fCutCharge)
3369       return kTRUE;
3370     else if (fCutCharge && fCharge * track->GetSign() > 0)
3371       return kTRUE;
3372   }
3373   return kFALSE;
3374 }
3375
3376
3377 //-----------------------------------------------------------------------
3378  // part added by Natasha
3379 Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
3380 {
3381   //pid selection for heavy nuclei
3382   Bool_t select=kFALSE;
3383
3384   //if (!track) continue;
3385
3386   if (!track->GetInnerParam())
3387     return kFALSE;    //break;
3388
3389   const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
3390
3391   Double_t ptotTPC = tpcTrack->GetP();
3392   Double_t sigTPC = track->GetTPCsignal();
3393   Double_t dEdxBBA = 0.;
3394   Double_t dSigma = 0.;
3395
3396   switch (fParticleID)
3397   {
3398     case AliPID::kDeuteron:
3399       //pid=10;
3400       dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
3401           4.60e+00,
3402           8.9684e+00,
3403           1.640e-05,
3404           2.35e+00,
3405           2.35e+00);
3406       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
3407
3408       if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) )  )
3409       {select=kTRUE;}
3410       break;
3411
3412     case AliPID::kTriton:
3413       //pid=11;
3414       select=kFALSE;
3415       break;
3416
3417     case AliPID::kHe3:
3418       //pid=12;
3419       // ----- Pass 2 -------
3420       dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
3421           1.74962,
3422           27.4992,
3423           4.00313e-15,
3424           2.42485,
3425           8.31768);
3426       dSigma = (sigTPC -  dEdxBBA)/dEdxBBA;
3427       if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
3428       {select=kTRUE;}
3429       break;
3430
3431     case AliPID::kAlpha:
3432       //pid=13;
3433       select=kFALSE;
3434       break;
3435
3436     default:
3437       return kFALSE;
3438   }
3439
3440   return select;
3441 }
3442 // end part added by Natasha
3443 //-----------------------------------------------------------------------
3444 Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliAODTrack* track)
3445 {
3446     // do a simple combined cut on the n sigma from tpc and tof
3447     // with information of the pid response object (needs pid response task)
3448     // stub, not implemented yet
3449     if(!track) return kFALSE;
3450     return kFALSE;
3451
3452 }
3453 //-----------------------------------------------------------------------------
3454 Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliESDtrack* track)
3455 {
3456     // do a simple combined cut on the n sigma from tpc and tof
3457     // with information of the pid response object (needs pid response task)
3458     // stub, not implemented yet
3459     if(!track) return kFALSE;
3460     return kFALSE;
3461
3462 }
3463
3464 //-----------------------------------------------------------------------
3465 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
3466  //set priors for the bayesian pid selection
3467   fCurrCentr = centrCur;
3468
3469   fBinLimitPID[0] = 0.300000;
3470   fBinLimitPID[1] = 0.400000;
3471   fBinLimitPID[2] = 0.500000;
3472   fBinLimitPID[3] = 0.600000;
3473   fBinLimitPID[4] = 0.700000;
3474   fBinLimitPID[5] = 0.800000;
3475   fBinLimitPID[6] = 0.900000;
3476   fBinLimitPID[7] = 1.000000;
3477   fBinLimitPID[8] = 1.200000;
3478   fBinLimitPID[9] = 1.400000;
3479   fBinLimitPID[10] = 1.600000;
3480   fBinLimitPID[11] = 1.800000;
3481   fBinLimitPID[12] = 2.000000;
3482   fBinLimitPID[13] = 2.200000;
3483   fBinLimitPID[14] = 2.400000;
3484   fBinLimitPID[15] = 2.600000;
3485   fBinLimitPID[16] = 2.800000;
3486   fBinLimitPID[17] = 3.000000;
3487
3488   // 0-10%
3489   if(centrCur < 10){
3490       fC[0][0] = 0.005;
3491       fC[0][1] = 0.005;
3492       fC[0][2] = 1.0000;
3493       fC[0][3] = 0.0010;
3494       fC[0][4] = 0.0010;
3495
3496       fC[1][0] = 0.005;
3497       fC[1][1] = 0.005;
3498       fC[1][2] = 1.0000;
3499       fC[1][3] = 0.0168;
3500       fC[1][4] = 0.0122;
3501
3502       fC[2][0] = 0.005;
3503       fC[2][1] = 0.005;
3504       fC[2][2] = 1.0000;
3505       fC[2][3] = 0.0272;
3506       fC[2][4] = 0.0070;
3507
3508       fC[3][0] = 0.005;
3509       fC[3][1] = 0.005;
3510       fC[3][2] = 1.0000;
3511       fC[3][3] = 0.0562;
3512       fC[3][4] = 0.0258;
3513
3514       fC[4][0] = 0.005;
3515       fC[4][1] = 0.005;
3516       fC[4][2] = 1.0000;
3517       fC[4][3] = 0.0861;
3518       fC[4][4] = 0.0496;
3519
3520       fC[5][0] = 0.005;
3521       fC[5][1] = 0.005;
3522       fC[5][2] = 1.0000;
3523       fC[5][3] = 0.1168;
3524       fC[5][4] = 0.0740;
3525
3526       fC[6][0] = 0.005;
3527       fC[6][1] = 0.005;
3528       fC[6][2] = 1.0000;
3529       fC[6][3] = 0.1476;
3530       fC[6][4] = 0.0998;
3531
3532       fC[7][0] = 0.005;
3533       fC[7][1] = 0.005;
3534       fC[7][2] = 1.0000;
3535       fC[7][3] = 0.1810;
3536       fC[7][4] = 0.1296;
3537
3538       fC[8][0] = 0.005;
3539       fC[8][1] = 0.005;
3540       fC[8][2] = 1.0000;
3541       fC[8][3] = 0.2240;
3542       fC[8][4] = 0.1827;
3543
3544       fC[9][0] = 0.005;
3545       fC[9][1] = 0.005;
3546       fC[9][2] = 1.0000;
3547       fC[9][3] = 0.2812;
3548       fC[9][4] = 0.2699;
3549
3550       fC[10][0] = 0.005;
3551       fC[10][1] = 0.005;
3552       fC[10][2] = 1.0000;
3553       fC[10][3] = 0.3328;
3554       fC[10][4] = 0.3714;
3555
3556       fC[11][0] = 0.005;
3557       fC[11][1] = 0.005;
3558       fC[11][2] = 1.0000;
3559       fC[11][3] = 0.3780;
3560       fC[11][4] = 0.4810;
3561
3562       fC[12][0] = 0.005;
3563       fC[12][1] = 0.005;
3564       fC[12][2] = 1.0000;
3565       fC[12][3] = 0.4125;
3566       fC[12][4] = 0.5771;
3567
3568       fC[13][0] = 0.005;
3569       fC[13][1] = 0.005;
3570       fC[13][2] = 1.0000;
3571       fC[13][3] = 0.4486;
3572       fC[13][4] = 0.6799;
3573
3574       fC[14][0] = 0.005;
3575       fC[14][1] = 0.005;
3576       fC[14][2] = 1.0000;
3577       fC[14][3] = 0.4840;
3578       fC[14][4] = 0.7668;
3579
3580       fC[15][0] = 0.005;
3581       fC[15][1] = 0.005;
3582       fC[15][2] = 1.0000;
3583       fC[15][3] = 0.4971;
3584       fC[15][4] = 0.8288;
3585
3586       fC[16][0] = 0.005;
3587       fC[16][1] = 0.005;
3588       fC[16][2] = 1.0000;
3589       fC[16][3] = 0.4956;
3590       fC[16][4] = 0.8653;
3591
3592       fC[17][0] = 0.005;
3593       fC[17][1] = 0.005;
3594       fC[17][2] = 1.0000;
3595       fC[17][3] = 0.5173;
3596       fC[17][4] = 0.9059;
3597   }
3598   // 10-20%
3599   else if(centrCur < 20){
3600      fC[0][0] = 0.005;
3601       fC[0][1] = 0.005;
3602       fC[0][2] = 1.0000;
3603       fC[0][3] = 0.0010;
3604       fC[0][4] = 0.0010;
3605
3606       fC[1][0] = 0.005;
3607       fC[1][1] = 0.005;
3608       fC[1][2] = 1.0000;
3609       fC[1][3] = 0.0132;
3610       fC[1][4] = 0.0088;
3611
3612       fC[2][0] = 0.005;
3613       fC[2][1] = 0.005;
3614       fC[2][2] = 1.0000;
3615       fC[2][3] = 0.0283;
3616       fC[2][4] = 0.0068;
3617
3618       fC[3][0] = 0.005;
3619       fC[3][1] = 0.005;
3620       fC[3][2] = 1.0000;
3621       fC[3][3] = 0.0577;
3622       fC[3][4] = 0.0279;
3623
3624       fC[4][0] = 0.005;
3625       fC[4][1] = 0.005;
3626       fC[4][2] = 1.0000;
3627       fC[4][3] = 0.0884;
3628       fC[4][4] = 0.0534;
3629
3630       fC[5][0] = 0.005;
3631       fC[5][1] = 0.005;
3632       fC[5][2] = 1.0000;
3633       fC[5][3] = 0.1179;
3634       fC[5][4] = 0.0794;
3635
3636       fC[6][0] = 0.005;
3637       fC[6][1] = 0.005;
3638       fC[6][2] = 1.0000;
3639       fC[6][3] = 0.1480;
3640       fC[6][4] = 0.1058;
3641
3642       fC[7][0] = 0.005;
3643       fC[7][1] = 0.005;
3644       fC[7][2] = 1.0000;
3645       fC[7][3] = 0.1807;
3646       fC[7][4] = 0.1366;
3647
3648       fC[8][0] = 0.005;
3649       fC[8][1] = 0.005;
3650       fC[8][2] = 1.0000;
3651       fC[8][3] = 0.2219;
3652       fC[8][4] = 0.1891;
3653
3654       fC[9][0] = 0.005;
3655       fC[9][1] = 0.005;
3656       fC[9][2] = 1.0000;
3657       fC[9][3] = 0.2804;
3658       fC[9][4] = 0.2730;
3659
3660       fC[10][0] = 0.005;
3661       fC[10][1] = 0.005;
3662       fC[10][2] = 1.0000;
3663       fC[10][3] = 0.3283;
3664       fC[10][4] = 0.3660;
3665
3666       fC[11][0] = 0.005;
3667       fC[11][1] = 0.005;
3668       fC[11][2] = 1.0000;
3669       fC[11][3] = 0.3710;
3670       fC[11][4] = 0.4647;
3671
3672       fC[12][0] = 0.005;
3673       fC[12][1] = 0.005;
3674       fC[12][2] = 1.0000;
3675       fC[12][3] = 0.4093;
3676       fC[12][4] = 0.5566;
3677
3678       fC[13][0] = 0.005;
3679       fC[13][1] = 0.005;
3680       fC[13][2] = 1.0000;
3681       fC[13][3] = 0.4302;
3682       fC[13][4] = 0.6410;
3683
3684       fC[14][0] = 0.005;
3685       fC[14][1] = 0.005;
3686       fC[14][2] = 1.0000;
3687       fC[14][3] = 0.4649;
3688       fC[14][4] = 0.7055;
3689
3690       fC[15][0] = 0.005;
3691       fC[15][1] = 0.005;
3692       fC[15][2] = 1.0000;
3693       fC[15][3] = 0.4523;
3694       fC[15][4] = 0.7440;
3695
3696       fC[16][0] = 0.005;
3697       fC[16][1] = 0.005;
3698       fC[16][2] = 1.0000;
3699       fC[16][3] = 0.4591;
3700       fC[16][4] = 0.7799;
3701
3702       fC[17][0] = 0.005;
3703       fC[17][1] = 0.005;
3704       fC[17][2] = 1.0000;
3705       fC[17][3] = 0.4804;
3706       fC[17][4] = 0.8218;
3707   }
3708   // 20-30%
3709   else if(centrCur < 30){
3710      fC[0][0] = 0.005;
3711       fC[0][1] = 0.005;
3712       fC[0][2] = 1.0000;
3713       fC[0][3] = 0.0010;
3714       fC[0][4] = 0.0010;
3715
3716       fC[1][0] = 0.005;
3717       fC[1][1] = 0.005;
3718       fC[1][2] = 1.0000;
3719       fC[1][3] = 0.0102;
3720       fC[1][4] = 0.0064;
3721
3722       fC[2][0] = 0.005;
3723       fC[2][1] = 0.005;
3724       fC[2][2] = 1.0000;
3725       fC[2][3] = 0.0292;
3726       fC[2][4] = 0.0066;
3727
3728       fC[3][0] = 0.005;
3729       fC[3][1] = 0.005;
3730       fC[3][2] = 1.0000;
3731       fC[3][3] = 0.0597;
3732       fC[3][4] = 0.0296;
3733
3734       fC[4][0] = 0.005;
3735       fC[4][1] = 0.005;
3736       fC[4][2] = 1.0000;
3737       fC[4][3] = 0.0900;
3738       fC[4][4] = 0.0589;
3739
3740       fC[5][0] = 0.005;
3741       fC[5][1] = 0.005;
3742       fC[5][2] = 1.0000;
3743       fC[5][3] = 0.1199;
3744       fC[5][4] = 0.0859;
3745
3746       fC[6][0] = 0.005;
3747       fC[6][1] = 0.005;
3748       fC[6][2] = 1.0000;
3749       fC[6][3] = 0.1505;
3750       fC[6][4] = 0.1141;
3751
3752       fC[7][0] = 0.005;
3753       fC[7][1] = 0.005;
3754       fC[7][2] = 1.0000;
3755       fC[7][3] = 0.1805;
3756       fC[7][4] = 0.1454;
3757
3758       fC[8][0] = 0.005;
3759       fC[8][1] = 0.005;
3760       fC[8][2] = 1.0000;
3761       fC[8][3] = 0.2221;
3762       fC[8][4] = 0.2004;
3763
3764       fC[9][0] = 0.005;
3765       fC[9][1] = 0.005;
3766       fC[9][2] = 1.0000;
3767       fC[9][3] = 0.2796;
3768       fC[9][4] = 0.2838;
3769
3770       fC[10][0] = 0.005;
3771       fC[10][1] = 0.005;
3772       fC[10][2] = 1.0000;
3773       fC[10][3] = 0.3271;
3774       fC[10][4] = 0.3682;
3775
3776       fC[11][0] = 0.005;
3777       fC[11][1] = 0.005;
3778       fC[11][2] = 1.0000;
3779       fC[11][3] = 0.3648;
3780       fC[11][4] = 0.4509;
3781
3782       fC[12][0] = 0.005;
3783       fC[12][1] = 0.005;
3784       fC[12][2] = 1.0000;
3785       fC[12][3] = 0.3988;
3786       fC[12][4] = 0.5339;
3787
3788       fC[13][0] = 0.005;
3789       fC[13][1] = 0.005;
3790       fC[13][2] = 1.0000;
3791       fC[13][3] = 0.4315;
3792       fC[13][4] = 0.5995;
3793
3794       fC[14][0] = 0.005;
3795       fC[14][1] = 0.005;
3796       fC[14][2] = 1.0000;
3797       fC[14][3] = 0.4548;
3798       fC[14][4] = 0.6612;
3799
3800       fC[15][0] = 0.005;
3801       fC[15][1] = 0.005;
3802       fC[15][2] = 1.0000;
3803       fC[15][3] = 0.4744;
3804       fC[15][4] = 0.7060;
3805
3806       fC[16][0] = 0.005;
3807       fC[16][1] = 0.005;
3808       fC[16][2] = 1.0000;
3809       fC[16][3] = 0.4899;
3810       fC[16][4] = 0.7388;
3811
3812       fC[17][0] = 0.005;
3813       fC[17][1] = 0.005;
3814       fC[17][2] = 1.0000;
3815       fC[17][3] = 0.4411;
3816       fC[17][4] = 0.7293;
3817   }
3818   // 30-40%
3819   else if(centrCur < 40){
3820       fC[0][0] = 0.005;
3821       fC[0][1] = 0.005;
3822       fC[0][2] = 1.0000;
3823       fC[0][3] = 0.0010;
3824       fC[0][4] = 0.0010;
3825
3826       fC[1][0] = 0.005;
3827       fC[1][1] = 0.005;
3828       fC[1][2] = 1.0000;
3829       fC[1][3] = 0.0102;
3830       fC[1][4] = 0.0048;
3831
3832       fC[2][0] = 0.005;
3833       fC[2][1] = 0.005;
3834       fC[2][2] = 1.0000;
3835       fC[2][3] = 0.0306;
3836       fC[2][4] = 0.0079;
3837
3838       fC[3][0] = 0.005;
3839       fC[3][1] = 0.005;
3840       fC[3][2] = 1.0000;
3841       fC[3][3] = 0.0617;
3842       fC[3][4] = 0.0338;
3843
3844       fC[4][0] = 0.005;
3845       fC[4][1] = 0.005;
3846       fC[4][2] = 1.0000;
3847       fC[4][3] = 0.0920;
3848       fC[4][4] = 0.0652;
3849
3850       fC[5][0] = 0.005;
3851       fC[5][1] = 0.005;
3852       fC[5][2] = 1.0000;
3853       fC[5][3] = 0.1211;
3854       fC[5][4] = 0.0955;
3855
3856       fC[6][0] = 0.005;
3857       fC[6][1] = 0.005;
3858       fC[6][2] = 1.0000;
3859       fC[6][3] = 0.1496;
3860       fC[6][4] = 0.1242;
3861
3862       fC[7][0] = 0.005;
3863       fC[7][1] = 0.005;
3864       fC[7][2] = 1.0000;
3865       fC[7][3] = 0.1807;
3866       fC[7][4] = 0.1576;
3867
3868       fC[8][0] = 0.005;
3869       fC[8][1] = 0.005;
3870       fC[8][2] = 1.0000;
3871       fC[8][3] = 0.2195;
3872       fC[8][4] = 0.2097;
3873
3874       fC[9][0] = 0.005;
3875       fC[9][1] = 0.005;
3876       fC[9][2] = 1.0000;
3877       fC[9][3] = 0.2732;
3878       fC[9][4] = 0.2884;
3879
3880       fC[10][0] = 0.005;
3881       fC[10][1] = 0.005;
3882       fC[10][2] = 1.0000;
3883       fC[10][3] = 0.3204;
3884       fC[10][4] = 0.3679;
3885
3886       fC[11][0] = 0.005;
3887       fC[11][1] = 0.005;
3888       fC[11][2] = 1.0000;
3889       fC[11][3] = 0.3564;
3890       fC[11][4] = 0.4449;
3891
3892       fC[12][0] = 0.005;
3893       fC[12][1] = 0.005;
3894       fC[12][2] = 1.0000;
3895       fC[12][3] = 0.3791;
3896       fC[12][4] = 0.5052;
3897
3898       fC[13][0] = 0.005;
3899       fC[13][1] = 0.005;
3900       fC[13][2] = 1.0000;
3901       fC[13][3] = 0.4062;
3902       fC[13][4] = 0.5647;
3903
3904       fC[14][0] = 0.005;
3905       fC[14][1] = 0.005;
3906       fC[14][2] = 1.0000;
3907       fC[14][3] = 0.4234;
3908       fC[14][4] = 0.6203;
3909
3910       fC[15][0] = 0.005;
3911       fC[15][1] = 0.005;
3912       fC[15][2] = 1.0000;
3913       fC[15][3] = 0.4441;
3914       fC[15][4] = 0.6381;
3915
3916       fC[16][0] = 0.005;
3917       fC[16][1] = 0.005;
3918       fC[16][2] = 1.0000;
3919       fC[16][3] = 0.4629;
3920       fC[16][4] = 0.6496;
3921
3922       fC[17][0] = 0.005;
3923       fC[17][1] = 0.005;
3924       fC[17][2] = 1.0000;
3925       fC[17][3] = 0.4293;
3926       fC[17][4] = 0.6491;
3927   }
3928   // 40-50%
3929   else if(centrCur < 50){
3930       fC[0][0] = 0.005;
3931       fC[0][1] = 0.005;
3932       fC[0][2] = 1.0000;
3933       fC[0][3] = 0.0010;
3934       fC[0][4] = 0.0010;
3935
3936       fC[1][0] = 0.005;
3937       fC[1][1] = 0.005;
3938       fC[1][2] = 1.0000;
3939       fC[1][3] = 0.0093;
3940       fC[1][4] = 0.0057;
3941
3942       fC[2][0] = 0.005;
3943       fC[2][1] = 0.005;
3944       fC[2][2] = 1.0000;
3945       fC[2][3] = 0.0319;
3946       fC[2][4] = 0.0075;
3947
3948       fC[3][0] = 0.005;
3949       fC[3][1] = 0.005;
3950       fC[3][2] = 1.0000;
3951       fC[3][3] = 0.0639;
3952       fC[3][4] = 0.0371;
3953
3954       fC[4][0] = 0.005;
3955       fC[4][1] = 0.005;
3956       fC[4][2] = 1.0000;
3957       fC[4][3] = 0.0939;
3958       fC[4][4] = 0.0725;
3959
3960       fC[5][0] = 0.005;
3961       fC[5][1] = 0.005;
3962       fC[5][2] = 1.0000;
3963       fC[5][3] = 0.1224;
3964       fC[5][4] = 0.1045;
3965
3966       fC[6][0] = 0.005;
3967       fC[6][1] = 0.005;
3968       fC[6][2] = 1.0000;
3969       fC[6][3] = 0.1520;
3970       fC[6][4] = 0.1387;
3971
3972       fC[7][0] = 0.005;
3973       fC[7][1] = 0.005;
3974       fC[7][2] = 1.0000;
3975       fC[7][3] = 0.1783;
3976       fC[7][4] = 0.1711;
3977
3978       fC[8][0] = 0.005;
3979       fC[8][1] = 0.005;
3980       fC[8][2] = 1.0000;
3981       fC[8][3] = 0.2202;
3982       fC[8][4] = 0.2269;
3983
3984       fC[9][0] = 0.005;
3985       fC[9][1] = 0.005;
3986       fC[9][2] = 1.0000;
3987       fC[9][3] = 0.2672;
3988       fC[9][4] = 0.2955;
3989
3990       fC[10][0] = 0.005;
3991       fC[10][1] = 0.005;
3992       fC[10][2] = 1.0000;
3993       fC[10][3] = 0.3191;
3994       fC[10][4] = 0.3676;
3995
3996       fC[11][0] = 0.005;
3997       fC[11][1] = 0.005;
3998       fC[11][2] = 1.0000;
3999       fC[11][3] = 0.3434;
4000       fC[11][4] = 0.4321;
4001
4002       fC[12][0] = 0.005;
4003       fC[12][1] = 0.005;
4004       fC[12][2] = 1.0000;
4005       fC[12][3] = 0.3692;
4006       fC[12][4] = 0.4879;
4007
4008       fC[13][0] = 0.005;
4009       fC[13][1] = 0.005;
4010       fC[13][2] = 1.0000;
4011       fC[13][3] = 0.3993;
4012       fC[13][4] = 0.5377;
4013
4014       fC[14][0] = 0.005;
4015       fC[14][1] = 0.005;
4016       fC[14][2] = 1.0000;
4017       fC[14][3] = 0.3818;
4018       fC[14][4] = 0.5547;
4019
4020       fC[15][0] = 0.005;
4021       fC[15][1] = 0.005;
4022       fC[15][2] = 1.0000;
4023       fC[15][3] = 0.4003;
4024       fC[15][4] = 0.5484;
4025
4026       fC[16][0] = 0.005;
4027       fC[16][1] = 0.005;
4028       fC[16][2] = 1.0000;
4029       fC[16][3] = 0.4281;
4030       fC[16][4] = 0.5383;
4031
4032       fC[17][0] = 0.005;
4033       fC[17][1] = 0.005;
4034       fC[17][2] = 1.0000;
4035       fC[17][3] = 0.3960;
4036       fC[17][4] = 0.5374;
4037   }
4038   // 50-60%
4039   else if(centrCur < 60){
4040       fC[0][0] = 0.005;
4041       fC[0][1] = 0.005;
4042       fC[0][2] = 1.0000;
4043       fC[0][3] = 0.0010;
4044       fC[0][4] = 0.0010;
4045
4046       fC[1][0] = 0.005;
4047       fC[1][1] = 0.005;
4048       fC[1][2] = 1.0000;
4049       fC[1][3] = 0.0076;
4050       fC[1][4] = 0.0032;
4051
4052       fC[2][0] = 0.005;
4053       fC[2][1] = 0.005;
4054       fC[2][2] = 1.0000;
4055       fC[2][3] = 0.0329;
4056       fC[2][4] = 0.0085;
4057
4058       fC[3][0] = 0.005;
4059       fC[3][1] = 0.005;
4060       fC[3][2] = 1.0000;
4061       fC[3][3] = 0.0653;
4062       fC[3][4] = 0.0423;
4063
4064       fC[4][0] = 0.005;
4065       fC[4][1] = 0.005;
4066       fC[4][2] = 1.0000;
4067       fC[4][3] = 0.0923;
4068       fC[4][4] = 0.0813;
4069
4070       fC[5][0] = 0.005;
4071       fC[5][1] = 0.005;
4072       fC[5][2] = 1.0000;
4073       fC[5][3] = 0.1219;
4074       fC[5][4] = 0.1161;
4075
4076       fC[6][0] = 0.005;
4077       fC[6][1] = 0.005;
4078       fC[6][2] = 1.0000;
4079       fC[6][3] = 0.1519;
4080       fC[6][4] = 0.1520;
4081
4082       fC[7][0] = 0.005;
4083       fC[7][1] = 0.005;
4084       fC[7][2] = 1.0000;
4085       fC[7][3] = 0.1763;
4086       fC[7][4] = 0.1858;
4087
4088       fC[8][0] = 0.005;
4089       fC[8][1] = 0.005;
4090       fC[8][2] = 1.0000;
4091       fC[8][3] = 0.2178;
4092       fC[8][4] = 0.2385;
4093
4094       fC[9][0] = 0.005;
4095       fC[9][1] = 0.005;
4096       fC[9][2] = 1.0000;
4097       fC[9][3] = 0.2618;
4098       fC[9][4] = 0.3070;
4099
4100       fC[10][0] = 0.005;
4101       fC[10][1] = 0.005;
4102       fC[10][2] = 1.0000;
4103       fC[10][3] = 0.3067;
4104       fC[10][4] = 0.3625;
4105
4106       fC[11][0] = 0.005;
4107       fC[11][1] = 0.005;
4108       fC[11][2] = 1.0000;
4109       fC[11][3] = 0.3336;
4110       fC[11][4] = 0.4188;
4111
4112       fC[12][0] = 0.005;
4113       fC[12][1] = 0.005;
4114       fC[12][2] = 1.0000;
4115       fC[12][3] = 0.3706;
4116       fC[12][4] = 0.4511;
4117
4118       fC[13][0] = 0.005;
4119       fC[13][1] = 0.005;
4120       fC[13][2] = 1.0000;
4121       fC[13][3] = 0.3765;
4122       fC[13][4] = 0.4729;
4123
4124       fC[14][0] = 0.005;
4125       fC[14][1] = 0.005;
4126       fC[14][2] = 1.0000;
4127       fC[14][3] = 0.3942;
4128       fC[14][4] = 0.4855;
4129
4130       fC[15][0] = 0.005;
4131       fC[15][1] = 0.005;
4132       fC[15][2] = 1.0000;
4133       fC[15][3] = 0.4051;
4134       fC[15][4] = 0.4762;
4135
4136       fC[16][0] = 0.005;
4137       fC[16][1] = 0.005;
4138       fC[16][2] = 1.0000;
4139       fC[16][3] = 0.3843;
4140       fC[16][4] = 0.4763;
4141
4142       fC[17][0] = 0.005;
4143       fC[17][1] = 0.005;
4144       fC[17][2] = 1.0000;
4145       fC[17][3] = 0.4237;
4146       fC[17][4] = 0.4773;
4147   }
4148   // 60-70%
4149   else if(centrCur < 70){
4150          fC[0][0] = 0.005;
4151       fC[0][1] = 0.005;
4152       fC[0][2] = 1.0000;
4153       fC[0][3] = 0.0010;
4154       fC[0][4] = 0.0010;
4155
4156       fC[1][0] = 0.005;
4157       fC[1][1] = 0.005;
4158       fC[1][2] = 1.0000;
4159       fC[1][3] = 0.0071;
4160       fC[1][4] = 0.0012;
4161
4162       fC[2][0] = 0.005;
4163       fC[2][1] = 0.005;
4164       fC[2][2] = 1.0000;
4165       fC[2][3] = 0.0336;
4166       fC[2][4] = 0.0097;
4167
4168       fC[3][0] = 0.005;
4169       fC[3][1] = 0.005;
4170       fC[3][2] = 1.0000;
4171       fC[3][3] = 0.0662;
4172       fC[3][4] = 0.0460;
4173
4174       fC[4][0] = 0.005;
4175       fC[4][1] = 0.005;
4176       fC[4][2] = 1.0000;
4177       fC[4][3] = 0.0954;
4178       fC[4][4] = 0.0902;
4179
4180       fC[5][0] = 0.005;
4181       fC[5][1] = 0.005;
4182       fC[5][2] = 1.0000;
4183       fC[5][3] = 0.1181;
4184       fC[5][4] = 0.1306;
4185
4186       fC[6][0] = 0.005;
4187       fC[6][1] = 0.005;
4188       fC[6][2] = 1.0000;
4189       fC[6][3] = 0.1481;
4190       fC[6][4] = 0.1662;
4191
4192       fC[7][0] = 0.005;
4193       fC[7][1] = 0.005;
4194       fC[7][2] = 1.0000;
4195       fC[7][3] = 0.1765;
4196       fC[7][4] = 0.1963;
4197
4198       fC[8][0] = 0.005;
4199       fC[8][1] = 0.005;
4200       fC[8][2] = 1.0000;
4201       fC[8][3] = 0.2155;
4202       fC[8][4] = 0.2433;
4203
4204       fC[9][0] = 0.005;
4205       fC[9][1] = 0.005;
4206       fC[9][2] = 1.0000;
4207       fC[9][3] = 0.2580;
4208       fC[9][4] = 0.3022;
4209
4210       fC[10][0] = 0.005;
4211       fC[10][1] = 0.005;
4212       fC[10][2] = 1.0000;
4213       fC[10][3] = 0.2872;
4214       fC[10][4] = 0.3481;
4215
4216       fC[11][0] = 0.005;
4217       fC[11][1] = 0.005;
4218       fC[11][2] = 1.0000;
4219       fC[11][3] = 0.3170;
4220       fC[11][4] = 0.3847;
4221
4222       fC[12][0] = 0.005;
4223       fC[12][1] = 0.005;
4224       fC[12][2] = 1.0000;
4225       fC[12][3] = 0.3454;
4226       fC[12][4] = 0.4258;
4227
4228       fC[13][0] = 0.005;
4229       fC[13][1] = 0.005;
4230       fC[13][2] = 1.0000;
4231       fC[13][3] = 0.3580;
4232       fC[13][4] = 0.4299;
4233
4234       fC[14][0] = 0.005;
4235       fC[14][1] = 0.005;
4236       fC[14][2] = 1.0000;
4237       fC[14][3] = 0.3903;
4238       fC[14][4] = 0.4326;
4239
4240       fC[15][0] = 0.005;
4241       fC[15][1] = 0.005;
4242       fC[15][2] = 1.0000;
4243       fC[15][3] = 0.3690;
4244       fC[15][4] = 0.4491;
4245
4246       fC[16][0] = 0.005;
4247       fC[16][1] = 0.005;
4248       fC[16][2] = 1.0000;
4249       fC[16][3] = 0.4716;
4250       fC[16][4] = 0.4298;
4251
4252       fC[17][0] = 0.005;
4253       fC[17][1] = 0.005;
4254       fC[17][2] = 1.0000;
4255       fC[17][3] = 0.3875;
4256       fC[17][4] = 0.4083;
4257   }
4258   // 70-80%
4259   else if(centrCur < 80){
4260       fC[0][0] = 0.005;
4261       fC[0][1] = 0.005;
4262       fC[0][2] = 1.0000;
4263       fC[0][3] = 0.0010;
4264       fC[0][4] = 0.0010;
4265
4266       fC[1][0] = 0.005;
4267       fC[1][1] = 0.005;
4268       fC[1][2] = 1.0000;
4269       fC[1][3] = 0.0075;
4270       fC[1][4] = 0.0007;
4271
4272       fC[2][0] = 0.005;
4273       fC[2][1] = 0.005;
4274       fC[2][2] = 1.0000;
4275       fC[2][3] = 0.0313;
4276       fC[2][4] = 0.0124;
4277
4278       fC[3][0] = 0.005;
4279       fC[3][1] = 0.005;
4280       fC[3][2] = 1.0000;
4281       fC[3][3] = 0.0640;
4282       fC[3][4] = 0.0539;
4283
4284       fC[4][0] = 0.005;
4285       fC[4][1] = 0.005;
4286       fC[4][2] = 1.0000;
4287       fC[4][3] = 0.0923;
4288       fC[4][4] = 0.0992;
4289
4290       fC[5][0] = 0.005;
4291       fC[5][1] = 0.005;
4292       fC[5][2] = 1.0000;
4293       fC[5][3] = 0.1202;
4294       fC[5][4] = 0.1417;
4295
4296       fC[6][0] = 0.005;
4297       fC[6][1] = 0.005;
4298       fC[6][2] = 1.0000;
4299       fC[6][3] = 0.1413;
4300       fC[6][4] = 0.1729;
4301
4302       fC[7][0] = 0.005;
4303       fC[7][1] = 0.005;
4304       fC[7][2] = 1.0000;
4305       fC[7][3] = 0.1705;
4306       fC[7][4] = 0.1999;
4307
4308       fC[8][0] = 0.005;
4309       fC[8][1] = 0.005;
4310       fC[8][2] = 1.0000;
4311       fC[8][3] = 0.2103;
4312       fC[8][4] = 0.2472;
4313
4314       fC[9][0] = 0.005;
4315       fC[9][1] = 0.005;
4316       fC[9][2] = 1.0000;
4317       fC[9][3] = 0.2373;
4318       fC[9][4] = 0.2916;
4319
4320       fC[10][0] = 0.005;
4321       fC[10][1] = 0.005;
4322       fC[10][2] = 1.0000;
4323       fC[10][3] = 0.2824;
4324       fC[10][4] = 0.3323;
4325
4326       fC[11][0] = 0.005;
4327       fC[11][1] = 0.005;
4328       fC[11][2] = 1.0000;
4329       fC[11][3] = 0.3046;
4330       fC[11][4] = 0.3576;
4331
4332       fC[12][0] = 0.005;
4333       fC[12][1] = 0.005;
4334       fC[12][2] = 1.0000;
4335       fC[12][3] = 0.3585;
4336       fC[12][4] = 0.4003;
4337
4338       fC[13][0] = 0.005;
4339       fC[13][1] = 0.005;
4340       fC[13][2] = 1.0000;
4341       fC[13][3] = 0.3461;
4342       fC[13][4] = 0.3982;
4343
4344       fC[14][0] = 0.005;
4345       fC[14][1] = 0.005;
4346       fC[14][2] = 1.0000;
4347       fC[14][3] = 0.3362;
4348       fC[14][4] = 0.3776;
4349
4350       fC[15][0] = 0.005;
4351       fC[15][1] = 0.005;
4352       fC[15][2] = 1.0000;
4353       fC[15][3] = 0.3071;
4354       fC[15][4] = 0.3500;
4355
4356       fC[16][0] = 0.005;
4357       fC[16][1] = 0.005;
4358       fC[16][2] = 1.0000;
4359       fC[16][3] = 0.2914;
4360       fC[16][4] = 0.3937;
4361
4362       fC[17][0] = 0.005;
4363       fC[17][1] = 0.005;
4364       fC[17][2] = 1.0000;
4365       fC[17][3] = 0.3727;
4366       fC[17][4] = 0.3877;
4367   }
4368   // 80-100%
4369   else{
4370       fC[0][0] = 0.005;
4371       fC[0][1] = 0.005;
4372       fC[0][2] = 1.0000;
4373       fC[0][3] = 0.0010;
4374       fC[0][4] = 0.0010;
4375
4376       fC[1][0] = 0.005;
4377       fC[1][1] = 0.005;
4378       fC[1][2] = 1.0000;
4379       fC[1][3] = 0.0060;
4380       fC[1][4] = 0.0035;
4381
4382       fC[2][0] = 0.005;
4383       fC[2][1] = 0.005;
4384       fC[2][2] = 1.0000;
4385       fC[2][3] = 0.0323;
4386       fC[2][4] = 0.0113;
4387
4388       fC[3][0] = 0.005;
4389       fC[3][1] = 0.005;
4390       fC[3][2] = 1.0000;
4391       fC[3][3] = 0.0609;
4392       fC[3][4] = 0.0653;
4393
4394       fC[4][0] = 0.005;
4395       fC[4][1] = 0.005;
4396       fC[4][2] = 1.0000;
4397       fC[4][3] = 0.0922;
4398       fC[4][4] = 0.1076;
4399
4400       fC[5][0] = 0.005;
4401       fC[5][1] = 0.005;
4402       fC[5][2] = 1.0000;
4403       fC[5][3] = 0.1096;
4404       fC[5][4] = 0.1328;
4405
4406       fC[6][0] = 0.005;
4407       fC[6][1] = 0.005;
4408       fC[6][2] = 1.0000;
4409       fC[6][3] = 0.1495;
4410       fC[6][4] = 0.1779;
4411
4412       fC[7][0] = 0.005;
4413       fC[7][1] = 0.005;
4414       fC[7][2] = 1.0000;
4415       fC[7][3] = 0.1519;
4416       fC[7][4] = 0.1989;
4417
4418       fC[8][0] = 0.005;
4419       fC[8][1] = 0.005;
4420       fC[8][2] = 1.0000;
4421       fC[8][3] = 0.1817;
4422       fC[8][4] = 0.2472;
4423
4424       fC[9][0] = 0.005;
4425       fC[9][1] = 0.005;
4426       fC[9][2] = 1.0000;
4427       fC[9][3] = 0.2429;
4428       fC[9][4] = 0.2684;
4429
4430       fC[10][0] = 0.005;
4431       fC[10][1] = 0.005;
4432       fC[10][2] = 1.0000;
4433       fC[10][3] = 0.2760;
4434       fC[10][4] = 0.3098;
4435
4436       fC[11][0] = 0.005;
4437       fC[11][1] = 0.005;
4438       fC[11][2] = 1.0000;
4439       fC[11][3] = 0.2673;
4440       fC[11][4] = 0.3198;
4441
4442       fC[12][0] = 0.005;
4443       fC[12][1] = 0.005;
4444       fC[12][2] = 1.0000;
4445       fC[12][3] = 0.3165;
4446       fC[12][4] = 0.3564;
4447
4448       fC[13][0] = 0.005;
4449       fC[13][1] = 0.005;
4450       fC[13][2] = 1.0000;
4451       fC[13][3] = 0.3526;
4452       fC[13][4] = 0.3011;
4453
4454       fC[14][0] = 0.005;
4455       fC[14][1] = 0.005;
4456       fC[14][2] = 1.0000;
4457       fC[14][3] = 0.3788;
4458       fC[14][4] = 0.3011;
4459
4460       fC[15][0] = 0.005;
4461       fC[15][1] = 0.005;
4462       fC[15][2] = 1.0000;
4463       fC[15][3] = 0.3788;
4464       fC[15][4] = 0.3011;
4465
4466       fC[16][0] = 0.005;
4467       fC[16][1] = 0.005;
4468       fC[16][2] = 1.0000;
4469       fC[16][3] = 0.3788;
4470       fC[16][4] = 0.3011;
4471
4472       fC[17][0] = 0.005;
4473       fC[17][1] = 0.005;
4474       fC[17][2] = 1.0000;
4475       fC[17][3] = 0.3788;
4476       fC[17][4] = 0.3011;
4477   }
4478
4479   for(Int_t i=18;i<fgkPIDptBin;i++){
4480     fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
4481     fC[i][0] = fC[17][0];
4482     fC[i][1] = fC[17][1];
4483     fC[i][2] = fC[17][2];
4484     fC[i][3] = fC[17][3];
4485     fC[i][4] = fC[17][4];
4486   }
4487 }
4488
4489 //---------------------------------------------------------------//
4490 Bool_t AliFlowTrackCuts::TPCTOFagree(const AliVTrack *track)
4491 {
4492   //check pid agreement between TPC and TOF
4493   Bool_t status = kFALSE;
4494
4495   const Float_t c = 2.99792457999999984e-02;
4496
4497   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
4498
4499
4500   Double_t exptimes[9];
4501   track->GetIntegratedTimes(exptimes);
4502
4503   Float_t dedx = track->GetTPCsignal();
4504
4505   Float_t p = track->P();
4506   Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
4507   Float_t tl = exptimes[0]*c; // to work both for ESD and AOD
4508
4509   Float_t betagammares =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4510
4511   Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
4512
4513 //  printf("betagamma1 = %f\n",betagamma1);
4514
4515   if(betagamma1 < 0.1) betagamma1 = 0.1;
4516
4517   if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
4518   else betagamma1 = 100;
4519
4520   Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
4521 //  printf("betagamma2 = %f\n",betagamma2);
4522
4523   if(betagamma2 < 0.1) betagamma2 = 0.1;
4524
4525   if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
4526   else betagamma2 = 100;
4527
4528
4529   Float_t momtpc=track->GetTPCmomentum();
4530
4531   for(Int_t i=0;i < 5;i++){
4532     Float_t resolutionTOF =  fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
4533     if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
4534       Float_t dedxExp = 0;
4535       if(i==0) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
4536       else if(i==1) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
4537       else if(i==2) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
4538       else if(i==3) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
4539       else if(i==4) dedxExp =  fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
4540
4541       Float_t resolutionTPC = 2;
4542       if(i==0) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
4543       else if(i==1) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
4544       else if(i==2) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
4545       else if(i==3) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
4546       else if(i==4) resolutionTPC =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4547
4548       if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
4549         status = kTRUE;
4550       }
4551     }
4552   }
4553
4554   Float_t bb1 =  fESDpid.GetTPCResponse().Bethe(betagamma1);
4555   Float_t bb2 =  fESDpid.GetTPCResponse().Bethe(betagamma2);
4556   Float_t bbM =  fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
4557
4558
4559   //  status = kFALSE;
4560   // for nuclei
4561   Float_t resolutionTOFpr =   fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4562   Float_t resolutionTPCpr =   fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4563   if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4564      status = kTRUE;
4565   }
4566   else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4567      status = kTRUE;
4568   }
4569   else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4570      status = kTRUE;
4571   }
4572
4573   return status;
4574 }
4575 // end part added by F. Noferini
4576 //-----------------------------------------------------------------------
4577
4578 //-----------------------------------------------------------------------
4579 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
4580 {
4581   //get the name of the particle id source
4582   switch (s)
4583   {
4584     case kTPCdedx:
4585       return "TPCdedx";
4586     case kTPCbayesian:
4587       return "TPCbayesian";
4588     case kTOFbeta:
4589       return "TOFbeta";
4590     case kTPCpid:
4591       return "TPCpid";
4592     case kTOFpid:
4593       return "TOFpid";
4594     case kTOFbayesian:
4595       return "TOFbayesianPID";
4596     case kTOFbetaSimple:
4597       return "TOFbetaSimple";
4598     case kTPCNuclei:
4599       return "TPCnuclei";
4600     case kTPCTOFNsigma:
4601       return "TPCTOFNsigma";
4602     default:
4603       return "NOPID";
4604   }
4605 }
4606
4607 //-----------------------------------------------------------------------
4608 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
4609 {
4610   //return the name of the selected parameter type
4611   switch (type)
4612   {
4613     case kMC:
4614       return "MC";
4615     case kGlobal:
4616       return "Global";
4617     case kTPCstandalone:
4618       return "TPCstandalone";
4619     case kSPDtracklet:
4620       return "SPDtracklets";
4621     case kPMD:
4622       return "PMD";
4623     case kVZERO:
4624       return "VZERO";
4625     case kMUON:       // XZhang 20120604
4626       return "MUON";  // XZhang 20120604
4627     case kKink:
4628       return "Kink";
4629     case kV0:
4630       return "V0";
4631     default:
4632       return "unknown";
4633   }
4634 }
4635
4636 //-----------------------------------------------------------------------
4637 Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
4638 {
4639   //check PMD specific cuts
4640   //clean up from last iteration, and init label
4641   Int_t   det   = track->GetDetector();
4642   //Int_t   smn   = track->GetSmn();
4643   Float_t clsX  = track->GetClusterX();
4644   Float_t clsY  = track->GetClusterY();
4645   Float_t clsZ  = track->GetClusterZ();
4646   Float_t ncell = track->GetClusterCells();
4647   Float_t adc   = track->GetClusterADC();
4648
4649   fTrack = NULL;
4650   fMCparticle=NULL;
4651   fTrackLabel=-996;
4652
4653   fTrackEta = GetPmdEta(clsX,clsY,clsZ);
4654   fTrackPhi = GetPmdPhi(clsX,clsY);
4655   fTrackWeight = 1.0;
4656
4657   Bool_t pass=kTRUE;
4658   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
4659   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
4660   if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
4661   if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
4662   if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
4663
4664   return pass;
4665 }
4666
4667 //-----------------------------------------------------------------------
4668 Bool_t AliFlowTrackCuts::PassesVZEROcuts(Int_t id)
4669 {
4670   //check VZERO cuts
4671   if (id<0) return kFALSE;
4672
4673   //clean up from last iter
4674   ClearTrack();
4675
4676   fTrackPhi = TMath::PiOver4()*(0.5+id%8);
4677
4678   // 10102013 weighting vzero tiles - rbertens@cern.ch
4679   if(!fVZEROgainEqualization) {
4680       // if for some reason the equalization is not initialized (e.g. 2011 data)
4681       // the fVZEROxpol[] weights are used to enable or disable vzero rings
4682     if(id<32) {   // v0c side
4683       fTrackEta = -3.45+0.5*(id/8);
4684       if(id < 8) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[0];
4685       else if (id < 16 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[1];
4686       else if (id < 24 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[2];
4687       else if (id < 32 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[3];
4688     } else {       // v0a side
4689       fTrackEta = +4.8-0.6*((id/8)-4);
4690       if( id < 40) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[0];
4691       else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[1];
4692       else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[2];
4693       else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[3];
4694     }
4695   } else { // the equalization is initialized
4696      // note that disabled rings have already been excluded on calibration level in
4697      // AliFlowEvent (so for a disabled ring, fVZEROxpol is zero
4698     if(id<32) {    // v0c side
4699       fTrackEta = -3.45+0.5*(id/8);
4700       if(id < 8) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4701       else if (id < 16 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4702       else if (id < 24 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4703       else if (id < 32 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[3]/fVZEROgainEqualization->GetBinContent(1+id);
4704     } else {       // v0a side
4705       fTrackEta = +4.8-0.6*((id/8)-4);
4706       if( id < 40) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4707       else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4708       else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4709       else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[3]/fVZEROgainEqualization->GetBinContent(1+id);
4710     }
4711     // printf ( " tile %i and weight %.2f \n", id, fTrackWeight);
4712   }
4713
4714   if (fLinearizeVZEROresponse && id < 64)
4715   {
4716     //this is only needed in pass1 of LHC10h
4717     Float_t multVZERO[fgkNumberOfVZEROtracks];
4718     Float_t dummy=0.0;
4719     AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multVZERO);
4720     fTrackWeight = multVZERO[id];
4721   }
4722
4723   Bool_t pass=kTRUE;
4724   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
4725   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
4726
4727   return pass;
4728 }
4729
4730 //-----------------------------------------------------------------------------
4731 Bool_t AliFlowTrackCuts::PassesMuonCuts(AliVParticle* vparticle)
4732 {
4733 // XZhang 20120604
4734   ClearTrack();
4735   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
4736   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
4737   else fMCparticle=NULL;
4738
4739   AliESDMuonTrack *esdTrack = dynamic_cast<AliESDMuonTrack*>(vparticle);
4740   AliAODTrack     *aodTrack = dynamic_cast<AliAODTrack*>(vparticle);
4741   if ((!esdTrack) && (!aodTrack)) return kFALSE;
4742   if (!fMuonTrackCuts->IsSelected(vparticle)) return kFALSE;
4743   HandleVParticle(vparticle);   if (!fTrack)  return kFALSE;
4744   return kTRUE;
4745 }
4746
4747 //----------------------------------------------------------------------------//
4748 Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
4749 {
4750   //get PMD "track" eta coordinate
4751   Float_t rpxpy, theta, eta;
4752   rpxpy  = TMath::Sqrt(xPos*xPos + yPos*yPos);
4753   theta  = TMath::ATan2(rpxpy,zPos);
4754   eta    = -TMath::Log(TMath::Tan(0.5*theta));
4755   return eta;
4756 }
4757
4758 //--------------------------------------------------------------------------//
4759 Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
4760 {
4761   //Get PMD "track" phi coordinate
4762   Float_t pybypx, phi = 0., phi1;
4763   if(xPos==0)
4764     {
4765       if(yPos>0) phi = 90.;
4766       if(yPos<0) phi = 270.;
4767     }
4768   if(xPos != 0)
4769     {
4770       pybypx = yPos/xPos;
4771       if(pybypx < 0) pybypx = - pybypx;
4772       phi1 = TMath::ATan(pybypx)*180./3.14159;
4773
4774       if(xPos > 0 && yPos > 0) phi = phi1;        // 1st Quadrant
4775       if(xPos < 0 && yPos > 0) phi = 180 - phi1;  // 2nd Quadrant
4776       if(xPos < 0 && yPos < 0) phi = 180 + phi1;  // 3rd Quadrant
4777       if(xPos > 0 && yPos < 0) phi = 360 - phi1;  // 4th Quadrant
4778
4779     }
4780   phi = phi*3.14159/180.;
4781   return   phi;
4782 }
4783
4784 //---------------------------------------------------------------//
4785 void AliFlowTrackCuts::Browse(TBrowser* b)
4786 {
4787   //some browsing capabilities
4788   if (fQA) b->Add(fQA);
4789 }
4790
4791 //---------------------------------------------------------------//
4792 Long64_t AliFlowTrackCuts::Merge(TCollection* list)
4793 {
4794   //merge
4795   if (!fQA || !list) return 0;
4796   if (list->IsEmpty()) return 0;
4797   AliFlowTrackCuts* obj=NULL;
4798   TList tmplist;
4799   TIter next(list);
4800   while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
4801   {
4802     if (obj==this) continue;
4803     tmplist.Add(obj->GetQA());
4804   }
4805   return fQA->Merge(&tmplist);
4806 }
4807
4808
4809