PWG2/FLOW/AliFlowTasks/AliFlowTrackCuts.cxx

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