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