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