PWG/FLOW/Tasks/AliFlowEventCuts.cxx

   1  /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   2  *                                                                        *
   3  * Author: The ALICE Off-line Project.                                    *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 /* $Id$ */
  16
  17 // AliFlowEventCuts:
  18 // An event cut class for the flow framework
  19 //
  20 // origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
  21
  22 #include <limits.h>
  23 #include <float.h>
  24 #include <TList.h>
  25 #include <TH1F.h>
  26 #include <TH2F.h>
  27 #include <TBrowser.h>
  28 #include "TMath.h"
  29 #include "TNamed.h"
  30 #include "AliVVertex.h"
  31 #include "AliVEvent.h"
  32 #include "AliESDEvent.h"
  33 #include "AliAODEvent.h"
  34 #include "AliAODHeader.h"
  35 #include "AliCentrality.h"
  36 #include "AliESDVZERO.h"
  37 #include "AliMultiplicity.h"
  38 #include "AliMCEvent.h"
  39 #include "AliFlowEventCuts.h"
  40 #include "AliFlowTrackCuts.h"
  41 #include "AliTriggerAnalysis.h"
  42 #include "AliCollisionGeometry.h"
  43 #include "AliGenEventHeader.h"
  44 #include <iostream>
  45 using namespace std;
  46
  47 ClassImp(AliFlowEventCuts)
  48
  49 //-----------------------------------------------------------------------
  50 AliFlowEventCuts::AliFlowEventCuts():
  51   AliFlowEventSimpleCuts(),
  52   fQA(NULL),
  53   fCutNumberOfTracks(kFALSE),
  54   fNumberOfTracksMax(INT_MAX),
  55   fNumberOfTracksMin(INT_MIN),
  56   fCutRefMult(kFALSE),
  57   fRefMultMethod(kTPConly),
  58   fUseAliESDtrackCutsRefMult(kFALSE),
  59   fRefMultMethodAliESDtrackCuts(AliESDtrackCuts::kTrackletsITSTPC),
  60   fRefMultMax(INT_MAX),
  61   fRefMultMin(INT_MIN),
  62   fRefMultCuts(NULL),
  63   fMeanPtCuts(NULL),
  64   fStandardTPCcuts(NULL),
  65   fStandardGlobalCuts(NULL),
  66   fCutPrimaryVertexX(kFALSE),
  67   fPrimaryVertexXmax(INT_MAX),
  68   fPrimaryVertexXmin(INT_MIN),
  69   fCutPrimaryVertexY(kFALSE),
  70   fPrimaryVertexYmax(INT_MAX),
  71   fPrimaryVertexYmin(INT_MIN),
  72   fCutPrimaryVertexZ(kFALSE),
  73   fPrimaryVertexZmax(INT_MAX),
  74   fPrimaryVertexZmin(INT_MIN),
  75   fCutNContributors(kFALSE),
  76   fNContributorsMax(INT_MAX),
  77   fNContributorsMin(INT_MIN),
  78   fCutMeanPt(kFALSE),
  79   fMeanPtMax(-DBL_MAX),
  80   fMeanPtMin(DBL_MAX),
  81   fCutSPDvertexerAnomaly(kFALSE),
  82   fCutSPDTRKVtxZ(kFALSE),
  83   fCutTPCmultiplicityOutliers(kFALSE),
  84   fCutTPCmultiplicityOutliersAOD(kFALSE),
  85   fUseCentralityUnchecked(kFALSE),
  86   fCentralityPercentileMethod(kTPConly),
  87   fCutZDCtiming(kFALSE),
  88   fTrigAna(),
  89   fCutImpactParameter(kFALSE),
  90   fImpactParameterMin(0.0),
  91   fImpactParameterMax(100.0),
  92   fhistTPCvsGlobalMult(0),
  93   fData2011(kFALSE)
  94 {
  95   //constructor
  96 }
  97
  98 //-----------------------------------------------------------------------
  99 AliFlowEventCuts::AliFlowEventCuts(const char* name, const char* title):
 100   AliFlowEventSimpleCuts(name, title),
 101   fQA(NULL),
 102   fCutNumberOfTracks(kFALSE),
 103   fNumberOfTracksMax(INT_MAX),
 104   fNumberOfTracksMin(INT_MIN),
 105   fCutRefMult(kFALSE),
 106   fRefMultMethod(kTPConly),
 107   fUseAliESDtrackCutsRefMult(kFALSE),
 108   fRefMultMethodAliESDtrackCuts(AliESDtrackCuts::kTrackletsITSTPC),
 109   fRefMultMax(INT_MAX),
 110   fRefMultMin(INT_MIN),
 111   fRefMultCuts(NULL),
 112   fMeanPtCuts(NULL),
 113   fStandardTPCcuts(AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()),
 114   fStandardGlobalCuts(AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()),
 115   fCutPrimaryVertexX(kFALSE),
 116   fPrimaryVertexXmax(INT_MAX),
 117   fPrimaryVertexXmin(INT_MIN),
 118   fCutPrimaryVertexY(kFALSE),
 119   fPrimaryVertexYmax(INT_MAX),
 120   fPrimaryVertexYmin(INT_MIN),
 121   fCutPrimaryVertexZ(kFALSE),
 122   fPrimaryVertexZmax(INT_MAX),
 123   fPrimaryVertexZmin(INT_MIN),
 124   fCutNContributors(kFALSE),
 125   fNContributorsMax(INT_MAX),
 126   fNContributorsMin(INT_MIN),
 127   fCutMeanPt(kFALSE),
 128   fMeanPtMax(-DBL_MAX),
 129   fMeanPtMin(DBL_MAX),
 130   fCutSPDvertexerAnomaly(kFALSE),
 131   fCutSPDTRKVtxZ(kFALSE),
 132   fCutTPCmultiplicityOutliers(kFALSE),
 133   fCutTPCmultiplicityOutliersAOD(kFALSE),
 134   fUseCentralityUnchecked(kFALSE),
 135   fCentralityPercentileMethod(kTPConly),
 136   fCutZDCtiming(kFALSE),
 137   fTrigAna(),
 138   fCutImpactParameter(kFALSE),
 139   fImpactParameterMin(0.0),
 140   fImpactParameterMax(100.0),
 141   fhistTPCvsGlobalMult(0),
 142   fData2011(kFALSE)
 143 {
 144   //constructor
 145 }
 146
 147 ////-----------------------------------------------------------------------
 148 AliFlowEventCuts::AliFlowEventCuts(const AliFlowEventCuts& that):
 149   AliFlowEventSimpleCuts(that),
 150   fQA(NULL),
 151   fCutNumberOfTracks(that.fCutNumberOfTracks),
 152   fNumberOfTracksMax(that.fNumberOfTracksMax),
 153   fNumberOfTracksMin(that.fNumberOfTracksMin),
 154   fCutRefMult(that.fCutRefMult),
 155   fRefMultMethod(that.fRefMultMethod),
 156   fUseAliESDtrackCutsRefMult(that.fUseAliESDtrackCutsRefMult),
 157   fRefMultMethodAliESDtrackCuts(that.fRefMultMethodAliESDtrackCuts),
 158   fRefMultMax(that.fRefMultMax),
 159   fRefMultMin(that.fRefMultMin),
 160   fRefMultCuts(NULL),
 161   fMeanPtCuts(NULL),
 162   fStandardTPCcuts(NULL),
 163   fStandardGlobalCuts(NULL),
 164   fCutPrimaryVertexX(that.fCutPrimaryVertexX),
 165   fPrimaryVertexXmax(that.fPrimaryVertexXmax),
 166   fPrimaryVertexXmin(that.fPrimaryVertexXmin),
 167   fCutPrimaryVertexY(that.fCutPrimaryVertexX),
 168   fPrimaryVertexYmax(that.fPrimaryVertexYmax),
 169   fPrimaryVertexYmin(that.fPrimaryVertexYmin),
 170   fCutPrimaryVertexZ(that.fCutPrimaryVertexX),
 171   fPrimaryVertexZmax(that.fPrimaryVertexZmax),
 172   fPrimaryVertexZmin(that.fPrimaryVertexZmin),
 173   fCutNContributors(that.fCutNContributors),
 174   fNContributorsMax(that.fNContributorsMax),
 175   fNContributorsMin(that.fNContributorsMin),
 176   fCutMeanPt(that.fCutMeanPt),
 177   fMeanPtMax(that.fMeanPtMax),
 178   fMeanPtMin(that.fMeanPtMin),
 179   fCutSPDvertexerAnomaly(that.fCutSPDvertexerAnomaly),
 180   fCutSPDTRKVtxZ(that.fCutSPDTRKVtxZ),
 181   fCutTPCmultiplicityOutliers(that.fCutTPCmultiplicityOutliers),
 182   fCutTPCmultiplicityOutliersAOD(that.fCutTPCmultiplicityOutliersAOD),
 183   fUseCentralityUnchecked(that.fUseCentralityUnchecked),
 184   fCentralityPercentileMethod(that.fCentralityPercentileMethod),
 185   fCutZDCtiming(that.fCutZDCtiming),
 186   fTrigAna(),
 187   fCutImpactParameter(that.fCutImpactParameter),
 188   fImpactParameterMin(that.fImpactParameterMin),
 189   fImpactParameterMax(that.fImpactParameterMax),
 190   fhistTPCvsGlobalMult(that.fhistTPCvsGlobalMult),
 191   fData2011(that.fData2011)
 192 {
 193   if (that.fQA) DefineHistograms();
 194   //copy constructor
 195   if (that.fRefMultCuts)
 196     fRefMultCuts = new AliFlowTrackCuts(*(that.fRefMultCuts));
 197   if (that.fMeanPtCuts)
 198     fMeanPtCuts = new AliFlowTrackCuts(*(that.fMeanPtCuts));
 199   fStandardTPCcuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
 200   fStandardGlobalCuts = AliFlowTrackCuts::GetStandardGlobalTrackCuts2010();
 201 }
 202
 203 ////-----------------------------------------------------------------------
 204 AliFlowEventCuts::~AliFlowEventCuts()
 205 {
 206   //dtor
 207   delete fMeanPtCuts;
 208   delete fRefMultCuts;
 209   delete fStandardGlobalCuts;
 210   delete fStandardTPCcuts;
 211   if (fQA) { fQA->SetOwner(); fQA->Delete(); delete fQA; }
 212 }
 213
 214 ////-----------------------------------------------------------------------
 215 AliFlowEventCuts& AliFlowEventCuts::operator=(const AliFlowEventCuts& that)
 216 {
 217   //assignment
 218   if (this==&that) return *this;
 219
 220   if (that.fQA)
 221   {
 222     if (fQA)
 223     {
 224       fQA->Delete();
 225       delete fQA;
 226     }
 227     fQA = static_cast<TList*>(that.fQA->Clone());
 228   }
 229   else
 230   {
 231     fQA->Delete();
 232     delete fQA;
 233     fQA=NULL;
 234   }
 235
 236   fCutNumberOfTracks=that.fCutNumberOfTracks;
 237   fNumberOfTracksMax=that.fNumberOfTracksMax;
 238   fNumberOfTracksMin=that.fNumberOfTracksMin;
 239   fCutRefMult=that.fCutRefMult;
 240   fRefMultMethod=that.fRefMultMethod;
 241   fUseAliESDtrackCutsRefMult=that.fUseAliESDtrackCutsRefMult;
 242   fRefMultMethodAliESDtrackCuts=that.fRefMultMethodAliESDtrackCuts;
 243   fRefMultMax=that.fRefMultMax;
 244   fRefMultMin=that.fRefMultMin;
 245   if (that.fRefMultCuts) *fRefMultCuts=*(that.fRefMultCuts);
 246   if (that.fMeanPtCuts) *fMeanPtCuts=*(that.fMeanPtCuts);
 247   fStandardTPCcuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010();
 248   fStandardGlobalCuts = AliFlowTrackCuts::GetStandardGlobalTrackCuts2010();
 249   fCutPrimaryVertexX=that.fCutPrimaryVertexX;
 250   fPrimaryVertexXmax=that.fPrimaryVertexXmax;
 251   fPrimaryVertexXmin=that.fPrimaryVertexXmin;
 252   fCutPrimaryVertexY=that.fCutPrimaryVertexY;
 253   fPrimaryVertexYmax=that.fPrimaryVertexYmax;
 254   fPrimaryVertexYmin=that.fPrimaryVertexYmin;
 255   fCutPrimaryVertexZ=that.fCutPrimaryVertexZ;
 256   fPrimaryVertexZmax=that.fPrimaryVertexZmax;
 257   fPrimaryVertexZmin=that.fPrimaryVertexZmin;
 258   fCutNContributors=that.fCutNContributors;
 259   fNContributorsMax=that.fNContributorsMax;
 260   fNContributorsMin=that.fNContributorsMin;
 261   fCutMeanPt=that.fCutMeanPt;
 262   fMeanPtMax=that.fMeanPtMax;
 263   fMeanPtMin=that.fMeanPtMin;
 264   fCutSPDvertexerAnomaly=that.fCutSPDvertexerAnomaly;
 265   fCutSPDTRKVtxZ=that.fCutSPDTRKVtxZ;
 266   fCutTPCmultiplicityOutliers=that.fCutTPCmultiplicityOutliers;
 267   fCutTPCmultiplicityOutliersAOD=that.fCutTPCmultiplicityOutliersAOD;
 268   fUseCentralityUnchecked=that.fUseCentralityUnchecked;
 269   fCentralityPercentileMethod=that.fCentralityPercentileMethod;
 270   fCutZDCtiming=that.fCutZDCtiming;
 271   fhistTPCvsGlobalMult=that.fhistTPCvsGlobalMult;
 272   fData2011=that.fData2011;
 273   return *this;
 274 }
 275
 276 //-----------------------------------------------------------------------
 277 Bool_t AliFlowEventCuts::IsSelected(TObject* obj, TObject* objmc)
 278 {
 279   //check cuts
 280   AliVEvent* vevent = dynamic_cast<AliVEvent*>(obj);
 281   AliMCEvent* mcevent = dynamic_cast<AliMCEvent*>(objmc);
 282   if (vevent) return PassesCuts(vevent,mcevent);;
 283   return kFALSE;  //when passed wrong type of object
 284 }
 285 //-----------------------------------------------------------------------
 286 Bool_t AliFlowEventCuts::PassesCuts(AliVEvent *event, AliMCEvent *mcevent)
 287 {
 288   ///check if event passes cuts
 289   const AliVVertex* pvtx=event->GetPrimaryVertex();
 290   Double_t pvtxx = pvtx->GetX();
 291   Double_t pvtxy = pvtx->GetY();
 292   Double_t pvtxz = pvtx->GetZ();
 293   Int_t ncontrib = pvtx->GetNContributors();
 294   Bool_t pass=kTRUE;
 295   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*>(event);
 296   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*>(event);
 297   Int_t multTPC = 0;
 298   Int_t multGlobal = 0;
 299   multTPC = fStandardTPCcuts->Count(event);
 300   multGlobal = fStandardGlobalCuts->Count(event);
 301   if (fQA)
 302   {
 303     QAbefore(0)->Fill(pvtxz);
 304     QAbefore(1)->Fill(multGlobal,multTPC);
 305   }
 306   if ( fCutTPCmultiplicityOutliers && esdevent )
 307   {
 308     //this is pretty slow as we check the event track by track twice
 309     //this cut will work for 2010 PbPb data and is dependent on
 310     //TPC and ITS reco efficiency (e.g. geometry, calibration etc)
 311     if(esdevent){
 312       if (multTPC > ( 23+1.216*multGlobal)) {pass=kFALSE;}
 313       if (multTPC < (-20+1.087*multGlobal)) {pass=kFALSE;}
 314     }
 315   }
 316
 317
 318   if(fCutTPCmultiplicityOutliersAOD && aodevent) {
 319       // the AliFlowTrackCuts::Count() function will not work here, since the correlation cut uses different
 320       // track cuts
 321       multTPC = 0; // tpc mult estimate
 322       Int_t multGlob = 0; // global multiplicity
 323       Int_t nGoodTracks(aodevent->GetNumberOfTracks());
 324       if(!fData2011) { // cut on outliers
 325           for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 326               AliAODTrack* trackAOD = aodevent->GetTrack(iTracks);
 327               if (!trackAOD) continue;
 328               if (!(trackAOD->TestFilterBit(1))) continue;
 329               if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
 330               multTPC++;
 331           }
 332           for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 333               AliAODTrack* trackAOD = aodevent->GetTrack(iTracks);
 334               if (!trackAOD) continue;
 335               if (!(trackAOD->TestFilterBit(16))) continue;
 336               if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
 337               Double_t b[2] = {-99., -99.};
 338               Double_t bCov[3] = {-99., -99., -99.};
 339               if (!(trackAOD->PropagateToDCA(aodevent->GetPrimaryVertex(), aodevent->GetMagneticField(), 100., b, bCov))) continue;
 340               if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 341               multGlob++;
 342           } //track loop
 343           if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
 344       }
 345       if(fData2011) { // cut on outliers
 346           for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 347               AliAODTrack* trackAOD = aodevent->GetTrack(iTracks);
 348               if (!trackAOD) continue;
 349               if (!(trackAOD->TestFilterBit(1))) continue;
 350               if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
 351               multTPC++;
 352           }
 353           for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 354               AliAODTrack* trackAOD = aodevent->GetTrack(iTracks);
 355               if (!trackAOD) continue;
 356               if (!(trackAOD->TestFilterBit(16))) continue;
 357               if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
 358               Double_t b[2] = {-99., -99.};
 359               Double_t bCov[3] = {-99., -99., -99.};
 360               if (!(trackAOD->PropagateToDCA(aodevent->GetPrimaryVertex(), aodevent->GetMagneticField(), 100., b, bCov))) continue;
 361               if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 362               multGlob++;
 363           } //track loop
 364           if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
 365       }
 366   }
 367
 368
 369   if (fCutNContributors)
 370   {
 371     if (ncontrib < fNContributorsMin || ncontrib >= fNContributorsMax) pass=kFALSE;
 372   }
 373   if (fCutPrimaryVertexX)
 374   {
 375     if (pvtxx < fPrimaryVertexXmin || pvtxx >= fPrimaryVertexXmax) pass=kFALSE;
 376   }
 377   if (fCutPrimaryVertexY)
 378   {
 379     if (pvtxy < fPrimaryVertexYmin || pvtxy >= fPrimaryVertexYmax) pass=kFALSE;
 380   }
 381   if (fCutPrimaryVertexZ)
 382   {
 383     if (pvtxz < fPrimaryVertexZmin || pvtxz >= fPrimaryVertexZmax)
 384       pass=kFALSE;
 385   }
 386   if (fCutCentralityPercentile&&esdevent)
 387   {
 388     AliCentrality* centr = esdevent->GetCentrality();
 389     if (fUseCentralityUnchecked)
 390     {
 391       if (!centr->IsEventInCentralityClassUnchecked( fCentralityPercentileMin,
 392                                                      fCentralityPercentileMax,
 393                                                      CentrMethName(fCentralityPercentileMethod) ))
 394       {
 395         pass=kFALSE;
 396       }
 397     }
 398     else
 399     {
 400       if (!centr->IsEventInCentralityClass( fCentralityPercentileMin,
 401                                             fCentralityPercentileMax,
 402                                             CentrMethName(fCentralityPercentileMethod) ))
 403       {
 404         pass=kFALSE;
 405       }
 406     }
 407   }
 408   if (fCutSPDvertexerAnomaly&&esdevent)
 409   {
 410     const AliESDVertex* sdpvertex = esdevent->GetPrimaryVertexSPD();
 411     if (sdpvertex->GetNContributors()<1) pass=kFALSE;
 412     if (sdpvertex->GetDispersion()>0.04) pass=kFALSE;
 413     if (sdpvertex->GetZRes()>0.25) pass=kFALSE;
 414     const AliESDVertex* tpcvertex = esdevent->GetPrimaryVertexTPC();
 415     if (tpcvertex->GetNContributors()<1) pass=kFALSE;
 416     const AliMultiplicity* tracklets = esdevent->GetMultiplicity();
 417     if (tpcvertex->GetNContributors()<(-10.0+0.25*tracklets->GetNumberOfITSClusters(0)))
 418     {
 419       pass=kFALSE;
 420     }
 421   }
 422   if (fCutZDCtiming&&esdevent)
 423   {
 424     if (!fTrigAna.ZDCTimeTrigger(esdevent))
 425     {
 426       pass=kFALSE;
 427     }
 428   }
 429   if(fCutNumberOfTracks) {if ( event->GetNumberOfTracks() < fNumberOfTracksMin ||
 430                                event->GetNumberOfTracks() >= fNumberOfTracksMax ) pass=kFALSE;}
 431   if((fCutRefMult&&mcevent)||(fCutRefMult&&esdevent))
 432   {
 433     //reference multiplicity still to be defined
 434     Double_t refMult = RefMult(event,mcevent);
 435     if (refMult < fRefMultMin || refMult >= fRefMultMax )
 436     {
 437       pass=kFALSE;
 438     }
 439   }
 440
 441   // Handles AOD event
 442   if(aodevent) {
 443     if(fCutSPDTRKVtxZ) {
 444       Double_t tVtxZ = aodevent->GetPrimaryVertex()->GetZ();
 445       Double_t tSPDVtxZ = aodevent->GetPrimaryVertexSPD()->GetZ();
 446       if( TMath::Abs(tVtxZ-tSPDVtxZ) > 0.5 ) pass = kFALSE;
 447     }
 448     AliCentrality* centr = aodevent->GetHeader()->GetCentralityP();
 449     if(fCutTPCmultiplicityOutliers){
 450       Double_t v0Centr  = centr->GetCentralityPercentile("V0M");
 451       Double_t trkCentr = centr->GetCentralityPercentile("TRK");
 452       if(TMath::Abs(v0Centr-trkCentr) > 5) pass = kFALSE;
 453     }
 454     if (fCutCentralityPercentile) {
 455       if (fUseCentralityUnchecked) {
 456         if (!centr->IsEventInCentralityClassUnchecked( fCentralityPercentileMin,
 457                                                        fCentralityPercentileMax,
 458                                                        CentrMethName(fCentralityPercentileMethod) )) {
 459           pass = kFALSE;
 460         }
 461       } else {
 462         if (!centr->IsEventInCentralityClass( fCentralityPercentileMin,
 463                                               fCentralityPercentileMax,
 464                                               CentrMethName(fCentralityPercentileMethod) )) {
 465           pass = kFALSE;
 466         }
 467       }
 468     }
 469   }
 470
 471   if (fCutMeanPt)
 472   {
 473     Float_t meanpt=0.0;
 474     Int_t ntracks=event->GetNumberOfTracks();
 475     Int_t nselected=0;
 476     for (Int_t i=0; i<ntracks; i++)
 477     {
 478       AliVParticle* track = event->GetTrack(i);
 479       if (!track) continue;
 480       Bool_t localpass=kTRUE;
 481       if (fMeanPtCuts) localpass=fMeanPtCuts->IsSelected(track);
 482       if (localpass)
 483       {
 484         meanpt += track->Pt();
 485         nselected++;
 486       }
 487     }
 488     if (nselected) meanpt=meanpt/nselected;
 489     if (meanpt<fMeanPtMin || meanpt >= fMeanPtMax) pass=kFALSE;
 490   }
 491
 492   //impact parameter cut
 493   if(fCutImpactParameter) {
 494     Double_t gImpactParameter = 0.;
 495     if(mcevent) {
 496       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(mcevent)->GenEventHeader());
 497       if(headerH)
 498         gImpactParameter = headerH->ImpactParameter();
 499     }
 500     if ((gImpactParameter < fImpactParameterMin) || (gImpactParameter >= fImpactParameterMax ))
 501       pass=kFALSE;
 502   }
 503
 504   if (fQA&&pass)
 505   {
 506     QAafter(1)->Fill(multGlobal,multTPC);
 507     QAafter(0)->Fill(pvtxz);
 508   }
 509   return pass;
 510 }
 511
 512 //-----------------------------------------------------------------------
 513 Float_t AliFlowEventCuts::GetCentrality(AliVEvent* event, AliMCEvent* /*mcEvent*/)
 514 {
 515   //get the centrality percentile of the event
 516   AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*>(event);
 517   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(event);
 518
 519   Float_t centrality=-1.;
 520
 521   AliCentrality* centr = NULL;
 522   if (esdEvent)
 523     centr = esdEvent->GetCentrality();
 524   if (aodEvent)
 525     centr = aodEvent->GetHeader()->GetCentralityP();
 526
 527   if (!centr) return -1.;
 528
 529   if (fUseCentralityUnchecked)
 530     centrality=centr->GetCentralityPercentileUnchecked(CentrMethName(fCentralityPercentileMethod));
 531   else
 532     centrality=centr->GetCentralityPercentile(CentrMethName(fCentralityPercentileMethod));
 533
 534   return centrality;
 535 }
 536
 537 //-----------------------------------------------------------------------
 538 const char* AliFlowEventCuts::CentrMethName(refMultMethod method) const
 539 {
 540   //get the string for refmultmethod, for use with AliCentrality in
 541   //the cut on centrality percentile
 542   switch (method)
 543   {
 544     case kSPDtracklets:
 545       return "TKL";
 546     case kSPD1clusters:
 547       return "CL1";
 548     case kTPConly:
 549       return "TRK";
 550     case kVZERO:
 551       return "V0M";
 552     default:
 553       return "";
 554   }
 555 }
 556 //-----------------------------------------------------------------------
 557 AliFlowEventCuts* AliFlowEventCuts::StandardCuts()
 558 {
 559   //make a set of standard event cuts, caller becomes owner
 560   AliFlowEventCuts* cuts = new AliFlowEventCuts();
 561   return cuts;
 562 }
 563
 564 //-----------------------------------------------------------------------
 565 Int_t AliFlowEventCuts::RefMult(AliVEvent* event, AliMCEvent *mcEvent)
 566 {
 567   //calculate the reference multiplicity, if all fails return 0
 568   AliESDVZERO* vzero = NULL;
 569   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*>(event);
 570
 571   if (fUseAliESDtrackCutsRefMult && esdevent)
 572   {
 573     //use the standard ALICE reference multiplicity with the default eta range
 574     return AliESDtrackCuts::GetReferenceMultiplicity(esdevent, fRefMultMethodAliESDtrackCuts);
 575   }
 576
 577   if (fRefMultMethod==kTPConly && !fRefMultCuts)
 578   {
 579     fRefMultCuts = AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts();
 580     fRefMultCuts->SetEtaRange(-0.8,0.8);
 581     fRefMultCuts->SetPtMin(0.15);
 582   }
 583   else if (fRefMultMethod==kSPDtracklets && !fRefMultCuts)
 584   {
 585     fRefMultCuts = new AliFlowTrackCuts("tracklet refmult cuts");
 586     fRefMultCuts->SetParamType(AliFlowTrackCuts::kSPDtracklet);
 587     fRefMultCuts->SetEtaRange(-0.8,0.8);
 588   }
 589   else if (fRefMultMethod==kVZERO)
 590   {
 591     if (!esdevent) return 0;
 592     vzero=esdevent->GetVZEROData();
 593     if (!vzero) return 0;
 594     return TMath::Nint(vzero->GetMTotV0A()+vzero->GetMTotV0C());
 595   }
 596   else if (fRefMultMethod==kSPD1clusters)
 597   {
 598     if (!esdevent) return 0;
 599     const AliMultiplicity* mult = esdevent->GetMultiplicity();
 600     if (!mult) return 0;
 601     return mult->GetNumberOfITSClusters(1);
 602   }
 603
 604   Int_t refmult=0;
 605   fRefMultCuts->SetEvent(event,mcEvent);
 606   Int_t numberOfInputObjects = fRefMultCuts->GetNumberOfInputObjects();
 607   for (Int_t i=0; i<numberOfInputObjects; i++) {
 608     if (fRefMultCuts->IsSelected(fRefMultCuts->GetInputObject(i),i))
 609       refmult++;
 610   }
 611   return refmult;
 612 }
 613 //_____________________________________________________________________________
 614 void AliFlowEventCuts::DefineHistograms()
 615 {
 616   //define QA histos
 617   if (fQA) return;
 618
 619   Bool_t adddirstatus = TH1::AddDirectoryStatus();
 620   TH1::AddDirectory(kFALSE);
 621   fQA = new TList(); fQA->SetOwner();
 622   fQA->SetName(Form("%s QA",GetName()));
 623   TList* before = new TList(); before->SetOwner();
 624   before->SetName("before");
 625   TList* after = new TList(); after->SetOwner();
 626   after->SetName("after");
 627   fQA->Add(before);
 628   fQA->Add(after);
 629   before->Add(new TH1F("zvertex",";z;event cout",500,-15.,15.)); //0
 630   after->Add(new TH1F("zvertex",";z;event cout",500,-15.,15.)); //0
 631   before->Add(new TH2F("fTPCvsGlobalMult","TPC only vs Global track multiplicity;global;TPC only",500,0,2500,500,0,3500));//1
 632   after->Add(new TH2F("fTPCvsGlobalMult","TPC only vs Global track multiplicity;global;TPC only",500,0,2500,500,0,3500));//1
 633   TH1::AddDirectory(adddirstatus);
 634 }
 635
 636 //---------------------------------------------------------------//
 637 void AliFlowEventCuts::Browse(TBrowser* b)
 638 {
 639   //some browsing capabilities
 640   if (fQA) b->Add(fQA);
 641 }
 642
 643 //---------------------------------------------------------------//
 644 Long64_t AliFlowEventCuts::Merge(TCollection* list)
 645 {
 646   //merge
 647   Int_t number=0;
 648   AliFlowEventCuts* obj;
 649   if (!list) return 0;
 650   if (list->GetEntries()<1) return 0;
 651   TIter next(list);
 652   while ( (obj = dynamic_cast<AliFlowEventCuts*>(next())) )
 653   {
 654     if (obj==this) continue;
 655     TList listwrapper;
 656     listwrapper.Add(obj->GetQA());
 657     fQA->Merge(&listwrapper);
 658     number++;
 659   }
 660   return number;
 661 }
 662
 663