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