PWGHF/vertexingHF/AliRDHFCuts.cxx

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