PWGHF/vertexingHF/AliRDHFCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, 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 /////////////////////////////////////////////////////////////
  19 //
  20 // Base class for cuts on AOD reconstructed heavy-flavour decay
  21 //
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24 #include <Riostream.h>
  25
  26 #include "AliVEvent.h"
  27 #include "AliESDEvent.h"
  28 #include "AliAODEvent.h"
  29 #include "AliVVertex.h"
  30 #include "AliESDVertex.h"
  31 #include "AliLog.h"
  32 #include "AliAODVertex.h"
  33 #include "AliESDtrack.h"
  34 #include "AliAODTrack.h"
  35 #include "AliESDtrackCuts.h"
  36 #include "AliCentrality.h"
  37 #include "AliAODRecoDecayHF.h"
  38 #include "AliAnalysisVertexingHF.h"
  39 #include "AliAODMCHeader.h"
  40 #include "AliAODMCParticle.h"
  41 #include "AliVertexerTracks.h"
  42 #include "AliRDHFCuts.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliInputEventHandler.h"
  45 #include "AliPIDResponse.h"
  46 #include "TRandom.h"
  47
  48 using std::cout;
  49 using std::endl;
  50
  51 ClassImp(AliRDHFCuts)
  52
  53 //--------------------------------------------------------------------------
  54 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
  55 AliAnalysisCuts(name,title),
  56 fMinVtxType(3),
  57 fMinVtxContr(1),
  58 fMaxVtxRedChi2(1e6),
  59 fMaxVtxZ(10.),
  60 fMinSPDMultiplicity(0),
  61 fTriggerMask(AliVEvent::kAnyINT),
  62 fUseOnlyOneTrigger(kFALSE),
  63 fTrackCuts(0),
  64 fnPtBins(1),
  65 fnPtBinLimits(1),
  66 fPtBinLimits(0),
  67 fnVars(1),
  68 fVarNames(0),
  69 fnVarsForOpt(0),
  70 fVarsForOpt(0),
  71 fGlobalIndex(1),
  72 fCutsRD(0),
  73 fIsUpperCut(0),
  74 fUsePID(kFALSE),
  75 fUseAOD049(kFALSE),
  76 fPidHF(0),
  77 fWhyRejection(0),
  78 fEvRejectionBits(0),
  79 fRemoveDaughtersFromPrimary(kFALSE),
  80 fRecomputePrimVertex(kFALSE),
  81 fUseMCVertex(kFALSE),
  82 fUsePhysicsSelection(kTRUE),
  83 fOptPileup(0),
  84 fMinContrPileup(3),
  85 fMinDzPileup(0.6),
  86 fUseCentrality(0),
  87 fMinCentrality(0.),
  88 fMaxCentrality(100.),
  89 fFixRefs(kFALSE),
  90 fIsSelectedCuts(0),
  91 fIsSelectedPID(0),
  92 fMinPtCand(-1.),
  93 fMaxPtCand(100000.),
  94 fKeepSignalMC(kFALSE),
  95 fIsCandTrackSPDFirst(kFALSE),
  96 fMaxPtCandTrackSPDFirst(0.),
  97 fApplySPDDeadPbPb2011(kFALSE),
  98 fRemoveTrackletOutliers(kFALSE),
  99 fCutOnzVertexSPD(0),
 100 fKinkReject(kFALSE),
 101   fUseTrackSelectionWithFilterBits(kTRUE),
 102   fUseCentrFlatteningInMC(kFALSE),
 103   fHistCentrDistr(0x0)
 104 {
 105   //
 106   // Default Constructor
 107   //
 108   fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
 109 }
 110 //--------------------------------------------------------------------------
 111 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
 112   AliAnalysisCuts(source),
 113   fMinVtxType(source.fMinVtxType),
 114   fMinVtxContr(source.fMinVtxContr),
 115   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
 116   fMaxVtxZ(source.fMaxVtxZ),
 117   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
 118   fTriggerMask(source.fTriggerMask),
 119   fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
 120   fTriggerClass(),
 121   fTrackCuts(0),
 122   fnPtBins(source.fnPtBins),
 123   fnPtBinLimits(source.fnPtBinLimits),
 124   fPtBinLimits(0),
 125   fnVars(source.fnVars),
 126   fVarNames(0),
 127   fnVarsForOpt(source.fnVarsForOpt),
 128   fVarsForOpt(0),
 129   fGlobalIndex(source.fGlobalIndex),
 130   fCutsRD(0),
 131   fIsUpperCut(0),
 132   fUsePID(source.fUsePID),
 133   fUseAOD049(source.fUseAOD049),
 134   fPidHF(0),
 135   fWhyRejection(source.fWhyRejection),
 136   fEvRejectionBits(source.fEvRejectionBits),
 137   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
 138   fRecomputePrimVertex(source.fRecomputePrimVertex),
 139   fUseMCVertex(source.fUseMCVertex),
 140   fUsePhysicsSelection(source.fUsePhysicsSelection),
 141   fOptPileup(source.fOptPileup),
 142   fMinContrPileup(source.fMinContrPileup),
 143   fMinDzPileup(source.fMinDzPileup),
 144   fUseCentrality(source.fUseCentrality),
 145   fMinCentrality(source.fMinCentrality),
 146   fMaxCentrality(source.fMaxCentrality),
 147   fFixRefs(source.fFixRefs),
 148   fIsSelectedCuts(source.fIsSelectedCuts),
 149   fIsSelectedPID(source.fIsSelectedPID),
 150   fMinPtCand(source.fMinPtCand),
 151   fMaxPtCand(source.fMaxPtCand),
 152   fKeepSignalMC(source.fKeepSignalMC),
 153   fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
 154   fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
 155   fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
 156   fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
 157   fCutOnzVertexSPD(source.fCutOnzVertexSPD),
 158   fKinkReject(source.fKinkReject),
 159   fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
 160   fUseCentrFlatteningInMC(source.fUseCentrFlatteningInMC),
 161   fHistCentrDistr(0x0)
 162 {
 163   //
 164   // Copy constructor
 165   //
 166   cout<<"Copy constructor"<<endl;
 167   fTriggerClass[0] = source.fTriggerClass[0];
 168   fTriggerClass[1] = source.fTriggerClass[1];
 169   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
 170   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 171   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 172   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 173   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 174   if(source.fPidHF) SetPidHF(source.fPidHF);
 175   if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
 176   PrintAll();
 177
 178 }
 179 //--------------------------------------------------------------------------
 180 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
 181 {
 182   //
 183   // assignment operator
 184   //
 185   if(&source == this) return *this;
 186
 187   AliAnalysisCuts::operator=(source);
 188
 189   fMinVtxType=source.fMinVtxType;
 190   fMinVtxContr=source.fMinVtxContr;
 191   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
 192   fMaxVtxZ=source.fMaxVtxZ;
 193   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
 194   fTriggerMask=source.fTriggerMask;
 195   fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
 196   fTriggerClass[0]=source.fTriggerClass[0];
 197   fTriggerClass[1]=source.fTriggerClass[1];
 198   fnPtBins=source.fnPtBins;
 199   fnPtBinLimits=source.fnPtBinLimits;
 200   fnVars=source.fnVars;
 201   fGlobalIndex=source.fGlobalIndex;
 202   fnVarsForOpt=source.fnVarsForOpt;
 203   fUsePID=source.fUsePID;
 204   fUseAOD049=source.fUseAOD049;
 205   if(fPidHF) delete fPidHF;
 206   fPidHF=new AliAODPidHF(*(source.GetPidHF()));
 207   fWhyRejection=source.fWhyRejection;
 208   fEvRejectionBits=source.fEvRejectionBits;
 209   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
 210   fRecomputePrimVertex=source.fRecomputePrimVertex;
 211   fUseMCVertex=source.fUseMCVertex;
 212   fUsePhysicsSelection=source.fUsePhysicsSelection;
 213   fOptPileup=source.fOptPileup;
 214   fMinContrPileup=source.fMinContrPileup;
 215   fMinDzPileup=source.fMinDzPileup;
 216   fUseCentrality=source.fUseCentrality;
 217   fMinCentrality=source.fMinCentrality;
 218   fMaxCentrality=source.fMaxCentrality;
 219   fFixRefs=source.fFixRefs;
 220   fIsSelectedCuts=source.fIsSelectedCuts;
 221   fIsSelectedPID=source.fIsSelectedPID;
 222   fMinPtCand=source.fMinPtCand;
 223   fMaxPtCand=source.fMaxPtCand;
 224   fKeepSignalMC=source.fKeepSignalMC;
 225   fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
 226   fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
 227   fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
 228   fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
 229   fCutOnzVertexSPD=source.fCutOnzVertexSPD;
 230   fKinkReject=source.fKinkReject;
 231   fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
 232   if(fHistCentrDistr) delete fHistCentrDistr;
 233   fUseCentrFlatteningInMC=source.fUseCentrFlatteningInMC;
 234   if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
 235
 236   if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
 237   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 238   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 239   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 240   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 241   PrintAll();
 242
 243   return *this;
 244 }
 245 //--------------------------------------------------------------------------
 246 AliRDHFCuts::~AliRDHFCuts() {
 247   //
 248   // Default Destructor
 249   //
 250   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
 251   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
 252   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
 253   if(fCutsRD) {
 254     delete [] fCutsRD;
 255     fCutsRD=0;
 256   }
 257   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
 258   if(fPidHF){
 259     delete fPidHF;
 260     fPidHF=0;
 261   }
 262   if(fHistCentrDistr)delete fHistCentrDistr;
 263 }
 264 //---------------------------------------------------------------------------
 265 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
 266   //
 267   // Centrality selection
 268   //
 269   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
 270     AliWarning("Centrality estimator not valid");
 271     return 3;
 272   }else{
 273     Float_t centvalue=GetCentrality((AliAODEvent*)event);
 274     if (centvalue<-998.){//-999 if no centralityP
 275       return 0;
 276     }else{
 277       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
 278         return 2;
 279       }
 280       // selection to flatten the centrality distribution
 281       if(fHistCentrDistr){
 282         if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
 283       }
 284     }
 285   }
 286   return 0;
 287 }
 288
 289
 290 //-------------------------------------------------
 291 void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
 292   // set the histo for centrality flattening
 293   // the centrality is flatten in the range minCentr,maxCentr
 294   // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
 295   //                positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
 296   //                negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
 297   // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
 298
 299   if(maxCentr<minCentr){
 300     AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
 301   }
 302
 303   if(fHistCentrDistr)delete fHistCentrDistr;
 304   fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
 305   fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
 306   Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
 307   Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
 308   fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
 309   Double_t ref=0.,bincont=0.,binrefwidth=1.;
 310   Int_t binref=0;
 311   if(TMath::Abs(centrRef)<0.0001){
 312     binref=fHistCentrDistr->GetMinimumBin();
 313     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
 314     ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
 315   }
 316   else if(centrRef>0.){
 317     binref=h->FindBin(centrRef);
 318     if(binref<1||binref>h->GetNbinsX()){
 319       AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
 320     }
 321     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
 322     ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
 323   }
 324   else{
 325     if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
 326     binref=fHistCentrDistr->GetMaximumBin();
 327     binrefwidth=fHistCentrDistr->GetBinWidth(binref);
 328     ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
 329   }
 330
 331   for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
 332     if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
 333       bincont=h->GetBinContent(j);
 334       fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
 335       fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
 336     }
 337     else{
 338       h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
 339     }
 340   }
 341
 342   fHistCentrDistr->SetBinContent(0,switchTRand);
 343   return;
 344
 345 }
 346
 347 //-------------------------------------------------
 348 Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
 349   //
 350   //  Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
 351   //  Can be faster if it was required that fHistCentrDistr covers
 352   //  exactly the desired centrality range (e.g. part of the lines below should be done during the
 353   // setting of the histo) and TH1::SetMinimum called
 354   //
 355
 356   if(!fHistCentrDistr) return kTRUE;
 357   // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
 358   //   if(maxbin>fHistCentrDistr->GetNbinsX()){
 359   //     AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
 360   //   }
 361
 362   Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
 363   Double_t bincont=fHistCentrDistr->GetBinContent(bin);
 364   Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
 365
 366   if(fHistCentrDistr->GetBinContent(0)<-0.9999){
 367     if(gRandom->Uniform(1.)<bincont)return kTRUE;
 368     return kFALSE;
 369   }
 370
 371   if(centDigits*100.<bincont)return kTRUE;
 372   return kFALSE;
 373
 374 }
 375 //---------------------------------------------------------------------------
 376 void AliRDHFCuts::SetupPID(AliVEvent *event) {
 377   // Set the PID response object in the AliAODPidHF
 378   // in case of old PID sets the TPC dE/dx BB parameterization
 379
 380   if(fPidHF){
 381     if(fPidHF->GetPidResponse()==0x0){
 382       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 383       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 384       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 385       fPidHF->SetPidResponse(pidResp);
 386     }
 387     if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
 388     if(fPidHF->GetOldPid()) {
 389
 390       Bool_t isMC=kFALSE;
 391       TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 392       if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
 393
 394       // pp, from LHC10d onwards
 395       if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
 396          event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
 397       // pp, 2011 low energy run
 398       if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
 399         fPidHF->SetppLowEn2011(kTRUE);
 400         fPidHF->SetOnePad(kFALSE);
 401       }
 402       // PbPb LHC10h
 403       if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
 404       // MC
 405       if(isMC) fPidHF->SetMC(kTRUE);
 406       if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
 407         fPidHF->SetMClowenpp2011(kTRUE);
 408       fPidHF->SetBetheBloch();
 409     }else{
 410       // check that AliPIDResponse object was properly set in case of using OADB
 411       if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
 412     }
 413   }
 414 }
 415 //---------------------------------------------------------------------------
 416 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
 417   //
 418   // Event selection
 419   //
 420   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
 421
 422
 423
 424   fWhyRejection=0;
 425   fEvRejectionBits=0;
 426   Bool_t accept=kTRUE;
 427
 428   if(fRecomputePrimVertex){
 429     Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
 430     if(!vertOK){
 431       fWhyRejection=6;
 432       return kFALSE;
 433     }
 434   }
 435
 436   // check if it's MC
 437   Bool_t isMC=kFALSE;
 438   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 439   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
 440
 441
 442   SetupPID(event);
 443
 444   // trigger class
 445   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
 446   // don't do for MC and for PbPb 2010 data
 447   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
 448     if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
 449        (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
 450       fWhyRejection=5;
 451       fEvRejectionBits+=1<<kNotSelTrigger;
 452       accept=kFALSE;
 453     }
 454   }
 455
 456   // TEMPORARY FIX FOR GetEvent
 457   Int_t nTracks=((AliAODEvent*)event)->GetNTracks();
 458   for(Int_t itr=0; itr<nTracks; itr++){
 459     AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
 460     tr->SetAODEvent((AliAODEvent*)event);
 461   }
 462
 463   // TEMPORARY FIX FOR REFERENCES
 464   // Fix references to daughter tracks
 465   //  if(fFixRefs) {
 466   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
 467   //    fixer->FixReferences((AliAODEvent*)event);
 468   //    delete fixer;
 469   //  }
 470   //
 471
 472
 473   // physics selection requirements
 474   if(fUsePhysicsSelection){
 475     Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
 476     if(!isSelected) {
 477       if(accept) fWhyRejection=7;
 478       fEvRejectionBits+=1<<kPhysicsSelection;
 479       accept=kFALSE;
 480     }else{
 481       if(fUseOnlyOneTrigger){
 482         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
 483           if(accept) fWhyRejection=7;
 484           fEvRejectionBits+=1<<kPhysicsSelection;
 485           accept=kFALSE;
 486         }
 487       }
 488     }
 489   }
 490
 491   // vertex requirements
 492
 493   const AliVVertex *vertex = event->GetPrimaryVertex();
 494
 495   if(!vertex){
 496     accept=kFALSE;
 497     fEvRejectionBits+=1<<kNoVertex;
 498   }else{
 499     TString title=vertex->GetTitle();
 500     if(title.Contains("Z") && fMinVtxType>1){
 501       accept=kFALSE;
 502       fEvRejectionBits+=1<<kNoVertex;
 503     }
 504     else if(title.Contains("3D") && fMinVtxType>2){
 505       accept=kFALSE;
 506       fEvRejectionBits+=1<<kNoVertex;
 507     }
 508     if(vertex->GetNContributors()<fMinVtxContr){
 509       accept=kFALSE;
 510       fEvRejectionBits+=1<<kTooFewVtxContrib;
 511     }
 512     if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
 513       fEvRejectionBits+=1<<kZVtxOutFid;
 514       if(accept) fWhyRejection=6;
 515       accept=kFALSE;
 516     }
 517   }
 518
 519   if(fCutOnzVertexSPD>0){
 520     const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
 521     if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
 522       accept=kFALSE;
 523       fEvRejectionBits+=1<<kBadSPDVertex;
 524     }else{
 525       if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
 526         fEvRejectionBits+=1<<kZVtxSPDOutFid;
 527         if(accept) fWhyRejection=6;
 528         accept=kFALSE;
 529       }
 530       if(fCutOnzVertexSPD==2 && vertex){
 531         if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
 532           fEvRejectionBits+=1<<kZVtxSPDOutFid;
 533           if(accept) fWhyRejection=6;
 534           accept=kFALSE;
 535         }
 536       }
 537     }
 538   }
 539
 540   // pile-up rejection
 541   if(fOptPileup==kRejectPileupEvent){
 542     Int_t cutc=(Int_t)fMinContrPileup;
 543     Double_t cutz=(Double_t)fMinDzPileup;
 544     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
 545       if(accept) fWhyRejection=1;
 546       fEvRejectionBits+=1<<kPileupSPD;
 547       accept=kFALSE;
 548     }
 549   }
 550
 551   // centrality selection
 552   if (fUseCentrality!=kCentOff) {
 553     Int_t rejection=IsEventSelectedInCentrality(event);
 554     Bool_t okCent=kFALSE;
 555     if(rejection==0) okCent=kTRUE;
 556     if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
 557     if(!okCent){
 558       if(accept) fWhyRejection=rejection;
 559       if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
 560       else fEvRejectionBits+=1<<kOutsideCentrality;
 561       accept=kFALSE;
 562     }
 563
 564   }
 565
 566   // PbPb2011 outliers in tracklets vs. VZERO
 567   if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
 568     if(fRemoveTrackletOutliers){
 569       Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
 570       Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
 571       Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
 572       if(ntracklets<1000. && v0cent<cutval){
 573         if(accept) fWhyRejection=2;
 574         fEvRejectionBits+=1<<kOutsideCentrality;
 575          accept=kFALSE;
 576       }
 577     }
 578   }
 579
 580   return accept;
 581 }
 582 //---------------------------------------------------------------------------
 583 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
 584   //
 585   // Daughter tracks selection
 586   //
 587   if(!fTrackCuts) return kTRUE;
 588
 589   Int_t ndaughters = d->GetNDaughters();
 590   AliAODVertex *vAOD = d->GetPrimaryVtx();
 591   Double_t pos[3],cov[6];
 592   vAOD->GetXYZ(pos);
 593   vAOD->GetCovarianceMatrix(cov);
 594   const AliESDVertex vESD(pos,cov,100.,100);
 595
 596   Bool_t retval=kTRUE;
 597
 598   for(Int_t idg=0; idg<ndaughters; idg++) {
 599     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
 600     if(!dgTrack) {retval = kFALSE; continue;}
 601     //printf("charge %d\n",dgTrack->Charge());
 602     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
 603
 604     if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
 605       { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
 606
 607     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
 608   }
 609
 610   return retval;
 611 }
 612 //---------------------------------------------------------------------------
 613 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
 614   //
 615   // Convert to ESDtrack, relate to vertex and check cuts
 616   //
 617   if(!cuts) return kTRUE;
 618
 619   if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
 620
 621   Bool_t retval=kTRUE;
 622
 623   // convert to ESD track here
 624   AliESDtrack esdTrack(track);
 625   // set the TPC cluster info
 626   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
 627   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
 628   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
 629   // needed to calculate the impact parameters
 630   esdTrack.RelateToVertex(primary,0.,3.);
 631
 632   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
 633
 634
 635   if(fKinkReject){
 636    AliAODVertex *maybeKink=track->GetProdVertex();
 637    if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
 638   }
 639
 640   if(fOptPileup==kRejectTracksFromPileupVertex){
 641     // to be implemented
 642     // we need either to have here the AOD Event,
 643     // or to have the pileup vertex object
 644   }
 645
 646   if(fApplySPDDeadPbPb2011){
 647     Bool_t deadSPDLay1PbPb2011[20][4]={
 648       {kTRUE,kTRUE,kTRUE,kTRUE},
 649       {kTRUE,kTRUE,kTRUE,kTRUE},
 650       {kTRUE,kTRUE,kTRUE,kTRUE},
 651       {kTRUE,kTRUE,kTRUE,kTRUE},
 652       {kTRUE,kTRUE,kTRUE,kTRUE},
 653       {kFALSE,kFALSE,kTRUE,kTRUE},
 654       {kTRUE,kTRUE,kFALSE,kFALSE},
 655       {kTRUE,kTRUE,kTRUE,kTRUE},
 656       {kFALSE,kFALSE,kTRUE,kTRUE},
 657       {kTRUE,kTRUE,kTRUE,kTRUE},
 658       {kTRUE,kTRUE,kFALSE,kFALSE},
 659       {kTRUE,kTRUE,kTRUE,kTRUE},
 660       {kFALSE,kFALSE,kFALSE,kFALSE},
 661       {kFALSE,kFALSE,kTRUE,kTRUE},
 662       {kFALSE,kFALSE,kFALSE,kFALSE},
 663       {kFALSE,kFALSE,kFALSE,kFALSE},
 664       {kTRUE,kTRUE,kTRUE,kTRUE},
 665       {kTRUE,kTRUE,kFALSE,kFALSE},
 666       {kFALSE,kFALSE,kFALSE,kFALSE},
 667       {kFALSE,kFALSE,kFALSE,kFALSE}
 668     };
 669     Bool_t deadSPDLay2PbPb2011[40][4]={
 670       {kTRUE,kTRUE,kTRUE,kTRUE},
 671       {kTRUE,kTRUE,kTRUE,kTRUE},
 672       {kTRUE,kTRUE,kTRUE,kTRUE},
 673       {kTRUE,kTRUE,kTRUE,kTRUE},
 674       {kTRUE,kTRUE,kTRUE,kTRUE},
 675       {kTRUE,kTRUE,kTRUE,kTRUE},
 676       {kTRUE,kTRUE,kTRUE,kTRUE},
 677       {kTRUE,kTRUE,kTRUE,kTRUE},
 678       {kTRUE,kTRUE,kTRUE,kTRUE},
 679       {kTRUE,kTRUE,kTRUE,kTRUE},
 680       {kTRUE,kTRUE,kTRUE,kTRUE},
 681       {kTRUE,kTRUE,kTRUE,kTRUE},
 682       {kFALSE,kFALSE,kFALSE,kFALSE},
 683       {kFALSE,kFALSE,kTRUE,kTRUE},
 684       {kTRUE,kTRUE,kTRUE,kTRUE},
 685       {kTRUE,kTRUE,kTRUE,kTRUE},
 686       {kTRUE,kTRUE,kFALSE,kFALSE},
 687       {kTRUE,kTRUE,kTRUE,kTRUE},
 688       {kTRUE,kTRUE,kTRUE,kTRUE},
 689       {kTRUE,kTRUE,kTRUE,kTRUE},
 690       {kFALSE,kFALSE,kFALSE,kFALSE},
 691       {kFALSE,kFALSE,kFALSE,kFALSE},
 692       {kTRUE,kTRUE,kTRUE,kTRUE},
 693       {kTRUE,kTRUE,kTRUE,kTRUE},
 694       {kFALSE,kFALSE,kFALSE,kFALSE},
 695       {kFALSE,kFALSE,kFALSE,kFALSE},
 696       {kTRUE,kTRUE,kTRUE,kTRUE},
 697       {kTRUE,kTRUE,kTRUE,kTRUE},
 698       {kFALSE,kFALSE,kFALSE,kFALSE},
 699       {kFALSE,kFALSE,kFALSE,kFALSE},
 700       {kFALSE,kFALSE,kFALSE,kFALSE},
 701       {kFALSE,kFALSE,kFALSE,kFALSE},
 702       {kTRUE,kTRUE,kTRUE,kTRUE},
 703       {kTRUE,kTRUE,kTRUE,kTRUE},
 704       {kTRUE,kTRUE,kFALSE,kFALSE},
 705       {kTRUE,kTRUE,kTRUE,kTRUE},
 706       {kFALSE,kFALSE,kFALSE,kFALSE},
 707       {kFALSE,kFALSE,kFALSE,kFALSE},
 708       {kFALSE,kFALSE,kFALSE,kFALSE},
 709       {kFALSE,kFALSE,kFALSE,kFALSE}
 710     };
 711     Double_t xyz1[3],xyz2[3];
 712     esdTrack.GetXYZAt(3.9,0.,xyz1);
 713     esdTrack.GetXYZAt(7.6,0.,xyz2);
 714     Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
 715     if(phi1<0) phi1+=2*TMath::Pi();
 716     Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
 717     Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
 718     if(phi2<0) phi2+=2*TMath::Pi();
 719     Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
 720     Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
 721     Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
 722     Bool_t lay1ok=kFALSE;
 723     if(mod1>=0 && mod1<4 && lad1<20){
 724       lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
 725     }
 726     Bool_t lay2ok=kFALSE;
 727     if(mod2>=0 && mod2<4 && lad2<40){
 728       lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
 729     }
 730     if(!lay1ok && !lay2ok) retval=kFALSE;
 731   }
 732
 733   return retval;
 734 }
 735 //---------------------------------------------------------------------------
 736 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
 737   // Set the pt bins
 738
 739   if(fPtBinLimits) {
 740     delete [] fPtBinLimits;
 741     fPtBinLimits = NULL;
 742     printf("Changing the pt bins\n");
 743   }
 744
 745   if(nPtBinLimits != fnPtBins+1){
 746     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
 747     SetNPtBins(nPtBinLimits-1);
 748   }
 749
 750   fnPtBinLimits = nPtBinLimits;
 751   SetGlobalIndex();
 752   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
 753   fPtBinLimits = new Float_t[fnPtBinLimits];
 754   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
 755
 756   return;
 757 }
 758 //---------------------------------------------------------------------------
 759 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
 760   // Set the variable names
 761
 762   if(fVarNames) {
 763     delete [] fVarNames;
 764     fVarNames = NULL;
 765     //printf("Changing the variable names\n");
 766   }
 767   if(nVars!=fnVars){
 768     printf("Wrong number of variables: it has to be %d\n",fnVars);
 769     return;
 770   }
 771   //fnVars=nVars;
 772   fVarNames = new TString[nVars];
 773   fIsUpperCut = new Bool_t[nVars];
 774   for(Int_t iv=0; iv<nVars; iv++) {
 775     fVarNames[iv] = varNames[iv];
 776     fIsUpperCut[iv] = isUpperCut[iv];
 777   }
 778
 779   return;
 780 }
 781 //---------------------------------------------------------------------------
 782 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
 783   // Set the variables to be used for cuts optimization
 784
 785   if(fVarsForOpt) {
 786     delete [] fVarsForOpt;
 787     fVarsForOpt = NULL;
 788     //printf("Changing the variables for cut optimization\n");
 789   }
 790
 791   if(nVars==0){//!=fnVars) {
 792     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
 793     return;
 794   }
 795
 796   fnVarsForOpt = 0;
 797   fVarsForOpt = new Bool_t[fnVars];
 798   for(Int_t iv=0; iv<fnVars; iv++) {
 799     fVarsForOpt[iv]=forOpt[iv];
 800     if(fVarsForOpt[iv]) fnVarsForOpt++;
 801   }
 802
 803   return;
 804 }
 805
 806 //---------------------------------------------------------------------------
 807 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
 808   //
 809   // set centrality estimator
 810   //
 811   fUseCentrality=flag;
 812   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
 813
 814   return;
 815 }
 816
 817
 818 //---------------------------------------------------------------------------
 819 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
 820   //
 821   // store the cuts
 822   //
 823   if(nVars!=fnVars) {
 824     printf("Wrong number of variables: it has to be %d\n",fnVars);
 825     AliFatal("exiting");
 826   }
 827   if(nPtBins!=fnPtBins) {
 828     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
 829     AliFatal("exiting");
 830   }
 831
 832   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 833
 834
 835   for(Int_t iv=0; iv<fnVars; iv++) {
 836
 837     for(Int_t ib=0; ib<fnPtBins; ib++) {
 838
 839       //check
 840       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
 841         cout<<"Overflow, exit..."<<endl;
 842         return;
 843       }
 844
 845       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
 846
 847     }
 848   }
 849   return;
 850 }
 851 //---------------------------------------------------------------------------
 852 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
 853   //
 854   // store the cuts
 855   //
 856   if(glIndex != fGlobalIndex){
 857     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
 858     AliFatal("exiting");
 859   }
 860   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 861
 862   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
 863     fCutsRD[iGl] = cutsRDGlob[iGl];
 864   }
 865   return;
 866 }
 867 //---------------------------------------------------------------------------
 868 void AliRDHFCuts::PrintAll() const {
 869   //
 870   // print all cuts values
 871   //
 872
 873   printf("Minimum vtx type %d\n",fMinVtxType);
 874   printf("Minimum vtx contr %d\n",fMinVtxContr);
 875   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
 876   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
 877   printf("Use PID %d  OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
 878   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
 879   printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
 880   printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
 881   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
 882   if(fOptPileup==1) printf(" -- Reject pileup event");
 883   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
 884   if(fUseCentrality>0) {
 885     TString estimator="";
 886     if(fUseCentrality==1) estimator = "V0";
 887     if(fUseCentrality==2) estimator = "Tracks";
 888     if(fUseCentrality==3) estimator = "Tracklets";
 889     if(fUseCentrality==4) estimator = "SPD clusters outer";
 890     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
 891   }
 892   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
 893
 894   if(fVarNames){
 895     cout<<"Array of variables"<<endl;
 896     for(Int_t iv=0;iv<fnVars;iv++){
 897       cout<<fVarNames[iv]<<"\t";
 898     }
 899     cout<<endl;
 900   }
 901   if(fVarsForOpt){
 902     cout<<"Array of optimization"<<endl;
 903     for(Int_t iv=0;iv<fnVars;iv++){
 904       cout<<fVarsForOpt[iv]<<"\t";
 905     }
 906     cout<<endl;
 907   }
 908   if(fIsUpperCut){
 909     cout<<"Array of upper/lower cut"<<endl;
 910    for(Int_t iv=0;iv<fnVars;iv++){
 911      cout<<fIsUpperCut[iv]<<"\t";
 912    }
 913    cout<<endl;
 914   }
 915   if(fPtBinLimits){
 916     cout<<"Array of ptbin limits"<<endl;
 917     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
 918       cout<<fPtBinLimits[ib]<<"\t";
 919     }
 920     cout<<endl;
 921   }
 922   if(fCutsRD){
 923     cout<<"Matrix of cuts"<<endl;
 924    for(Int_t iv=0;iv<fnVars;iv++){
 925      for(Int_t ib=0;ib<fnPtBins;ib++){
 926        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
 927      }
 928      cout<<endl;
 929    }
 930    cout<<endl;
 931   }
 932   return;
 933 }
 934
 935 //--------------------------------------------------------------------------
 936 void AliRDHFCuts::PrintTrigger() const{
 937   // print the trigger selection
 938
 939   printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
 940
 941   cout<<" Trigger selection pattern: ";
 942   if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
 943   if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
 944   if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
 945   if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
 946   if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
 947   if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
 948   if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
 949   if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
 950   if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
 951   if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
 952   cout << endl<< endl;
 953
 954 }
 955
 956 //---------------------------------------------------------------------------
 957 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
 958   //
 959   // get the cuts
 960   //
 961
 962   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
 963
 964
 965   Int_t iv,ib;
 966   if(!cutsRD) {
 967     //cout<<"Initialization..."<<endl;
 968     cutsRD=new Float_t*[fnVars];
 969     for(iv=0; iv<fnVars; iv++) {
 970       cutsRD[iv] = new Float_t[fnPtBins];
 971     }
 972   }
 973
 974   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
 975     GetVarPtIndex(iGlobal,iv,ib);
 976     cutsRD[iv][ib] = fCutsRD[iGlobal];
 977   }
 978
 979   return;
 980 }
 981
 982 //---------------------------------------------------------------------------
 983 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
 984   //
 985   // give the global index from variable and pt bin
 986   //
 987   return iPtBin*fnVars+iVar;
 988 }
 989
 990 //---------------------------------------------------------------------------
 991 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
 992   //
 993   //give the index of the variable and of the pt bin from the global index
 994   //
 995   iPtBin=(Int_t)iGlob/fnVars;
 996   iVar=iGlob%fnVars;
 997
 998   return;
 999 }
1000
1001 //---------------------------------------------------------------------------
1002 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1003   //
1004   //give the pt bin where the pt lies.
1005   //
1006   Int_t ptbin=-1;
1007   if(pt<fPtBinLimits[0])return ptbin;
1008   for (Int_t i=0;i<fnPtBins;i++){
1009     if(pt<fPtBinLimits[i+1]) {
1010       ptbin=i;
1011       break;
1012     }
1013   }
1014   return ptbin;
1015 }
1016 //-------------------------------------------------------------------
1017 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1018   //
1019   // Give the value of cut set for the variable iVar and the pt bin iPtBin
1020   //
1021   if(!fCutsRD){
1022     cout<<"Cuts not iniziaisez yet"<<endl;
1023     return 0;
1024   }
1025   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1026 }
1027 //-------------------------------------------------------------------
1028 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1029   //
1030   // Get centrality percentile
1031   //
1032
1033   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1034   if(mcArray) {fUseAOD049=kFALSE;}
1035
1036   AliAODHeader *header=aodEvent->GetHeader();
1037   AliCentrality *centrality=header->GetCentralityP();
1038   Float_t cent=-999.;
1039   Bool_t isSelRun=kFALSE;
1040   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1041   if(!centrality) return cent;
1042   else{
1043     if (estimator==kCentV0M){
1044       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1045       if(cent<0){
1046         Int_t quality = centrality->GetQuality();
1047         if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1048           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1049         }else{
1050           Int_t runnum=aodEvent->GetRunNumber();
1051           for(Int_t ir=0;ir<5;ir++){
1052             if(runnum==selRun[ir]){
1053               isSelRun=kTRUE;
1054               break;
1055             }
1056           }
1057           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1058         }
1059       }
1060
1061       //temporary fix for AOD049 outliers
1062       if(fUseAOD049&&cent>=0){
1063         Float_t v0=0;
1064         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1065         v0+=aodV0->GetMTotV0A();
1066         v0+=aodV0->GetMTotV0C();
1067         if(cent==0&&v0<19500)return -1;//filtering issue
1068         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1069         Float_t val= 1.30552 +  0.147931 * v0;
1070         Float_t tklSigma[101]={176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86, 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654, 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334, 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224, 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255, 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398, 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235, 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504, 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1071         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1072       }
1073     }
1074     else {
1075       if (estimator==kCentTRK) {
1076         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1077         if(cent<0){
1078           Int_t quality = centrality->GetQuality();
1079           if(quality<=1){
1080             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1081           }else{
1082             Int_t runnum=aodEvent->GetRunNumber();
1083             for(Int_t ir=0;ir<5;ir++){
1084               if(runnum==selRun[ir]){
1085                 isSelRun=kTRUE;
1086                 break;
1087               }
1088             }
1089             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1090           }
1091         }
1092       }
1093       else{
1094         if (estimator==kCentTKL){
1095           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1096           if(cent<0){
1097             Int_t quality = centrality->GetQuality();
1098             if(quality<=1){
1099               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1100             }else{
1101               Int_t runnum=aodEvent->GetRunNumber();
1102               for(Int_t ir=0;ir<5;ir++){
1103                 if(runnum==selRun[ir]){
1104                   isSelRun=kTRUE;
1105                   break;
1106             }
1107               }
1108               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1109             }
1110           }
1111         }
1112         else{
1113           if (estimator==kCentCL1){
1114             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1115             if(cent<0){
1116               Int_t quality = centrality->GetQuality();
1117               if(quality<=1){
1118                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1119               }else{
1120                 Int_t runnum=aodEvent->GetRunNumber();
1121                 for(Int_t ir=0;ir<5;ir++){
1122                   if(runnum==selRun[ir]){
1123                     isSelRun=kTRUE;
1124                     break;
1125                   }
1126                 }
1127                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1128               }
1129             }
1130           }
1131           else {
1132             AliWarning("Centrality estimator not valid");
1133
1134           }
1135         }
1136       }
1137     }
1138   }
1139   return cent;
1140 }
1141 //-------------------------------------------------------------------
1142 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1143   //
1144   // Compare two cuts objects
1145   //
1146
1147   Bool_t areEqual=kTRUE;
1148
1149   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1150
1151   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1152
1153   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1154
1155   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1156
1157   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1158
1159   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1160   if(fTrackCuts){
1161     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1162
1163     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1164
1165     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1166
1167     if(fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)!=obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)) {printf("ClusterReq SPD %d  %d\n",fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD),obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)); areEqual=kFALSE;}
1168   }
1169
1170   if(fCutsRD) {
1171    for(Int_t iv=0;iv<fnVars;iv++) {
1172      for(Int_t ib=0;ib<fnPtBins;ib++) {
1173        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1174          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1175          areEqual=kFALSE;
1176        }
1177      }
1178    }
1179   }
1180
1181   return areEqual;
1182 }
1183 //---------------------------------------------------------------------------
1184 void AliRDHFCuts::MakeTable() const {
1185   //
1186   // print cuts values in table format
1187   //
1188
1189         TString ptString = "pT range";
1190         if(fVarNames && fPtBinLimits && fCutsRD){
1191                 TString firstLine(Form("*       %-15s",ptString.Data()));
1192                 for (Int_t ivar=0; ivar<fnVars; ivar++){
1193                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
1194                         if (ivar == fnVars){
1195                                 firstLine+="*\n";
1196                         }
1197                 }
1198                 Printf("%s",firstLine.Data());
1199
1200                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1201                         TString line;
1202                         if (ipt==fnPtBins-1){
1203                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
1204                         }
1205                         else{
1206                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1207                         }
1208                         for (Int_t ivar=0; ivar<fnVars; ivar++){
1209                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1210                         }
1211                         Printf("%s",line.Data());
1212                 }
1213
1214         }
1215
1216
1217   return;
1218 }
1219 //--------------------------------------------------------------------------
1220 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1221                                         AliAODEvent *aod) const
1222 {
1223   //
1224   // Recalculate primary vertex without daughters
1225   //
1226
1227   if(!aod) {
1228     AliError("Can not remove daughters from vertex without AOD event");
1229     return 0;
1230   }
1231
1232   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1233   if(!recvtx){
1234     AliDebug(2,"Removal of daughter tracks failed");
1235     return kFALSE;
1236   }
1237
1238
1239   //set recalculed primary vertex
1240   d->SetOwnPrimaryVtx(recvtx);
1241   delete recvtx;
1242
1243   return kTRUE;
1244 }
1245 //--------------------------------------------------------------------------
1246 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1247 {
1248   //
1249   // Recalculate primary vertex without daughters
1250   //
1251
1252   if(!aod) {
1253     AliError("Can not get MC vertex without AOD event");
1254     return kFALSE;
1255   }
1256
1257   // load MC header
1258   AliAODMCHeader *mcHeader =
1259     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1260   if(!mcHeader) {
1261     AliError("Can not get MC vertex without AODMCHeader event");
1262     return kFALSE;
1263   }
1264   Double_t pos[3];
1265   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1266   mcHeader->GetVertex(pos);
1267   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1268
1269   if(!recvtx){
1270     AliDebug(2,"Removal of daughter tracks failed");
1271     return kFALSE;
1272   }
1273
1274   //set recalculed primary vertex
1275   d->SetOwnPrimaryVtx(recvtx);
1276
1277   d->RecalculateImpPars(recvtx,aod);
1278
1279   delete recvtx;
1280
1281   return kTRUE;
1282 }
1283 //--------------------------------------------------------------------------
1284 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1285                                      AliAODEvent *aod,
1286                                      AliAODVertex *origownvtx) const
1287 {
1288   //
1289   // Clean-up own primary vertex if needed
1290   //
1291
1292   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1293     d->UnsetOwnPrimaryVtx();
1294     if(origownvtx) {
1295       d->SetOwnPrimaryVtx(origownvtx);
1296       delete origownvtx; origownvtx=NULL;
1297     }
1298     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1299   } else {
1300     if(origownvtx) {
1301       delete origownvtx; origownvtx=NULL;
1302     }
1303   }
1304   return;
1305 }
1306 //--------------------------------------------------------------------------
1307 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1308 {
1309   //
1310   // Checks if this candidate is matched to MC signal
1311   //
1312
1313   if(!aod) return kFALSE;
1314
1315   // get the MC array
1316   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1317
1318   if(!mcArray) return kFALSE;
1319
1320   // try to match
1321   Int_t label = d->MatchToMC(pdg,mcArray);
1322
1323   if(label>=0) {
1324     //printf("MATCH!\n");
1325     return kTRUE;
1326   }
1327
1328   return kFALSE;
1329 }
1330
1331
1332 //--------------------------------------------------------------------------
1333 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1334   // recompute event primary vertex from AOD tracks
1335
1336    AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1337    vertexer->SetITSMode();
1338    vertexer->SetMinClusters(3);
1339
1340    AliAODVertex* pvtx=event->GetPrimaryVertex();
1341    if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1342      Float_t diamondcovxy[3];
1343      event->GetDiamondCovXY(diamondcovxy);
1344      Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1345      Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1346      AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1347      vertexer->SetVtxStart(diamond);
1348      delete diamond; diamond=NULL;
1349    }
1350
1351    AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1352    if(!vertexESD) return kFALSE;
1353    if(vertexESD->GetNContributors()<=0) {
1354      //AliDebug(2,"vertexing failed");
1355      delete vertexESD; vertexESD=NULL;
1356      return kFALSE;
1357    }
1358    delete vertexer; vertexer=NULL;
1359
1360    // convert to AliAODVertex
1361    Double_t pos[3],cov[6],chi2perNDF;
1362    vertexESD->GetXYZ(pos); // position
1363    vertexESD->GetCovMatrix(cov); //covariance matrix
1364    chi2perNDF = vertexESD->GetChi2toNDF();
1365    delete vertexESD; vertexESD=NULL;
1366
1367    pvtx->SetPosition(pos[0],pos[1],pos[2]);
1368    pvtx->SetChi2perNDF(chi2perNDF);
1369    pvtx->SetCovMatrix(cov);
1370
1371    return kTRUE;
1372 }