PWG3/vertexingHF/AliRDHFCuts.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 /* $Id$ */\r
  17 \r
  18 /////////////////////////////////////////////////////////////\r
  19 //\r
  20 // Base class for cuts on AOD reconstructed heavy-flavour decay\r
  21 //\r
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it\r
  23 /////////////////////////////////////////////////////////////\r
  24 #include <Riostream.h>\r
  25 \r
  26 #include "AliVEvent.h"\r
  27 #include "AliESDEvent.h"\r
  28 #include "AliAODEvent.h"\r
  29 #include "AliVVertex.h"\r
  30 #include "AliESDVertex.h"\r
  31 #include "AliLog.h"\r
  32 #include "AliAODVertex.h"\r
  33 #include "AliESDtrack.h"\r
  34 #include "AliAODTrack.h"\r
  35 #include "AliESDtrackCuts.h"\r
  36 #include "AliCentrality.h"\r
  37 #include "AliAODRecoDecayHF.h"\r
  38 #include "AliAnalysisVertexingHF.h"\r
  39 #include "AliAODMCHeader.h"\r
  40 #include "AliAODMCParticle.h"\r
  41 #include "AliRDHFCuts.h"\r
  42 #include "AliAnalysisManager.h"\r
  43 #include "AliInputEventHandler.h"\r
  44 #include "AliPIDResponse.h"\r
  45 \r
  46 ClassImp(AliRDHFCuts)\r
  47 \r
  48 //--------------------------------------------------------------------------\r
  49 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : \r
  50 AliAnalysisCuts(name,title),\r
  51 fMinVtxType(3),\r
  52 fMinVtxContr(1),\r
  53 fMaxVtxRedChi2(1e6),\r
  54 fMaxVtxZ(10.),\r
  55 fMinSPDMultiplicity(0),\r
  56 fTriggerMask(AliVEvent::kAnyINT),\r
  57 fUseOnlyOneTrigger(kFALSE),\r
  58 fTriggerClass("CINT1"),\r
  59 fTrackCuts(0),\r
  60 fnPtBins(1),\r
  61 fnPtBinLimits(1),\r
  62 fPtBinLimits(0),\r
  63 fnVars(1),\r
  64 fVarNames(0),\r
  65 fnVarsForOpt(0),\r
  66 fVarsForOpt(0),\r
  67 fGlobalIndex(1),\r
  68 fCutsRD(0),\r
  69 fIsUpperCut(0),\r
  70 fUsePID(kFALSE),\r
  71 fUseAOD049(kFALSE),\r
  72 fPidHF(0),\r
  73 fWhyRejection(0),\r
  74 fEvRejectionBits(0),\r
  75 fRemoveDaughtersFromPrimary(kFALSE),\r
  76 fUseMCVertex(kFALSE),\r
  77 fUsePhysicsSelection(kTRUE),\r
  78 fOptPileup(0),\r
  79 fMinContrPileup(3),\r
  80 fMinDzPileup(0.6),\r
  81 fUseCentrality(0),\r
  82 fMinCentrality(0.),\r
  83 fMaxCentrality(100.),\r
  84 fFixRefs(kFALSE),\r
  85 fIsSelectedCuts(0),\r
  86 fIsSelectedPID(0),\r
  87 fMinPtCand(-1.),\r
  88 fMaxPtCand(100000.),\r
  89 fKeepSignalMC(kFALSE),\r
  90 fIsCandTrackSPDFirst(kFALSE),\r
  91 fMaxPtCandTrackSPDFirst(0.)\r
  92 {\r
  93   //\r
  94   // Default Constructor\r
  95   //\r
  96 }\r
  97 //--------------------------------------------------------------------------\r
  98 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :\r
  99   AliAnalysisCuts(source),\r
 100   fMinVtxType(source.fMinVtxType),\r
 101   fMinVtxContr(source.fMinVtxContr),\r
 102   fMaxVtxRedChi2(source.fMaxVtxRedChi2),\r
 103   fMaxVtxZ(source.fMaxVtxZ),\r
 104   fMinSPDMultiplicity(source.fMinSPDMultiplicity),\r
 105   fTriggerMask(source.fTriggerMask),\r
 106   fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),\r
 107   fTriggerClass(source.fTriggerClass),\r
 108   fTrackCuts(0),\r
 109   fnPtBins(source.fnPtBins),\r
 110   fnPtBinLimits(source.fnPtBinLimits),\r
 111   fPtBinLimits(0),\r
 112   fnVars(source.fnVars),\r
 113   fVarNames(0),\r
 114   fnVarsForOpt(source.fnVarsForOpt),\r
 115   fVarsForOpt(0),\r
 116   fGlobalIndex(source.fGlobalIndex),\r
 117   fCutsRD(0),\r
 118   fIsUpperCut(0),\r
 119   fUsePID(source.fUsePID),\r
 120   fUseAOD049(source.fUseAOD049),\r
 121   fPidHF(0),\r
 122   fWhyRejection(source.fWhyRejection),\r
 123   fEvRejectionBits(source.fEvRejectionBits),\r
 124   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),\r
 125   fUseMCVertex(source.fUseMCVertex),\r
 126   fUsePhysicsSelection(source.fUsePhysicsSelection),\r
 127   fOptPileup(source.fOptPileup),\r
 128   fMinContrPileup(source.fMinContrPileup),\r
 129   fMinDzPileup(source.fMinDzPileup),\r
 130   fUseCentrality(source.fUseCentrality),\r
 131   fMinCentrality(source.fMinCentrality),\r
 132   fMaxCentrality(source.fMaxCentrality),\r
 133   fFixRefs(source.fFixRefs),\r
 134   fIsSelectedCuts(source.fIsSelectedCuts),\r
 135   fIsSelectedPID(source.fIsSelectedPID),\r
 136   fMinPtCand(source.fMinPtCand),\r
 137   fMaxPtCand(source.fMaxPtCand),\r
 138   fKeepSignalMC(source.fKeepSignalMC),\r
 139   fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),\r
 140   fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst)\r
 141 {\r
 142   //\r
 143   // Copy constructor\r
 144   //\r
 145   cout<<"Copy constructor"<<endl;\r
 146   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());\r
 147   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);\r
 148   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);\r
 149   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);\r
 150   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);\r
 151   if(source.fPidHF) SetPidHF(source.fPidHF);\r
 152   PrintAll();\r
 153 \r
 154 }\r
 155 //--------------------------------------------------------------------------\r
 156 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)\r
 157 {\r
 158   //\r
 159   // assignment operator\r
 160   //\r
 161   if(&source == this) return *this;\r
 162 \r
 163   AliAnalysisCuts::operator=(source);\r
 164 \r
 165   fMinVtxType=source.fMinVtxType;\r
 166   fMinVtxContr=source.fMinVtxContr;\r
 167   fMaxVtxRedChi2=source.fMaxVtxRedChi2;\r
 168   fMaxVtxZ=source.fMaxVtxZ;\r
 169   fMinSPDMultiplicity=source.fMinSPDMultiplicity;\r
 170   fTriggerMask=source.fTriggerMask;\r
 171   fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;\r
 172   fTriggerClass=source.fTriggerClass;\r
 173   fnPtBins=source.fnPtBins;\r
 174   fnPtBinLimits=source.fnPtBinLimits;\r
 175   fnVars=source.fnVars;\r
 176   fGlobalIndex=source.fGlobalIndex;\r
 177   fnVarsForOpt=source.fnVarsForOpt;\r
 178   fUsePID=source.fUsePID;\r
 179   fUseAOD049=source.fUseAOD049;\r
 180   SetPidHF(source.GetPidHF());\r
 181   fWhyRejection=source.fWhyRejection;\r
 182   fEvRejectionBits=source.fEvRejectionBits;\r
 183   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;\r
 184   fUseMCVertex=source.fUseMCVertex;\r
 185   fUsePhysicsSelection=source.fUsePhysicsSelection;\r
 186   fOptPileup=source.fOptPileup;\r
 187   fMinContrPileup=source.fMinContrPileup;\r
 188   fMinDzPileup=source.fMinDzPileup;\r
 189   fUseCentrality=source.fUseCentrality;\r
 190   fMinCentrality=source.fMinCentrality;\r
 191   fMaxCentrality=source.fMaxCentrality;\r
 192   fFixRefs=source.fFixRefs;\r
 193   fIsSelectedCuts=source.fIsSelectedCuts;\r
 194   fIsSelectedPID=source.fIsSelectedPID;\r
 195   fMinPtCand=source.fMinPtCand;\r
 196   fMaxPtCand=source.fMaxPtCand;\r
 197   fKeepSignalMC=source.fKeepSignalMC;\r
 198   fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;\r
 199   fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;\r
 200 \r
 201   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());\r
 202   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);\r
 203   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);\r
 204   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);\r
 205   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);\r
 206   PrintAll();\r
 207 \r
 208   return *this;\r
 209 }\r
 210 //--------------------------------------------------------------------------\r
 211 AliRDHFCuts::~AliRDHFCuts() {\r
 212   //  \r
 213   // Default Destructor\r
 214   //\r
 215   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}\r
 216   if(fVarNames) {delete [] fVarNames; fVarNames=0;}\r
 217   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}\r
 218   if(fCutsRD) {\r
 219     delete [] fCutsRD;\r
 220     fCutsRD=0;\r
 221   }\r
 222   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}\r
 223   if(fPidHF){ \r
 224     delete fPidHF;\r
 225     fPidHF=0;\r
 226   }\r
 227 }\r
 228 //---------------------------------------------------------------------------\r
 229 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {\r
 230   //\r
 231   // Centrality selection\r
 232   //\r
 233   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){    \r
 234     AliWarning("Centrality estimator not valid");    \r
 235     return 3;  \r
 236   }else{    \r
 237     Float_t centvalue=GetCentrality((AliAODEvent*)event);          \r
 238     if (centvalue<-998.){//-999 if no centralityP\r
 239       return 0;    \r
 240     }else{      \r
 241       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){\r
 242         return 2;      \r
 243       }    \r
 244     } \r
 245   }  \r
 246   return 0;\r
 247 }\r
 248 //---------------------------------------------------------------------------\r
 249 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {\r
 250   //\r
 251   // Event selection\r
 252   // \r
 253   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;\r
 254 \r
 255   fWhyRejection=0;\r
 256   fEvRejectionBits=0;\r
 257   Bool_t accept=kTRUE;\r
 258 \r
 259   // check if it's MC\r
 260   Bool_t isMC=kFALSE;\r
 261   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 262   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}\r
 263 \r
 264   // settings for the TPC dE/dx BB parameterization\r
 265   if(fPidHF && fPidHF->GetOldPid()) {\r
 266     // pp, from LHC10d onwards\r
 267     if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||\r
 268        event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);\r
 269     // pp, 2011 low energy run\r
 270     if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){\r
 271       fPidHF->SetppLowEn2011(kTRUE);\r
 272       fPidHF->SetOnePad(kFALSE);\r
 273     }\r
 274     // PbPb LHC10h\r
 275     if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);\r
 276     // MC\r
 277     if(isMC) fPidHF->SetMC(kTRUE);\r
 278     if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))\r
 279       fPidHF->SetMClowenpp2011(kTRUE);\r
 280   } \r
 281   else if(fPidHF && !fPidHF->GetOldPid()) {\r
 282    if(fPidHF->GetPidResponse()==0x0){\r
 283     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
 284     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();\r
 285     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();\r
 286     fPidHF->SetPidResponse(pidResp);\r
 287    }\r
 288   }\r
 289 \r
 290 \r
 291   // trigger class\r
 292   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();\r
 293   // don't do for MC and for PbPb 2010 data\r
 294   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {\r
 295     if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {\r
 296       fWhyRejection=5;\r
 297       fEvRejectionBits+=1<<kNotSelTrigger;\r
 298       accept=kFALSE;\r
 299     }\r
 300   }\r
 301 \r
 302   // TEMPORARY FIX FOR REFERENCES\r
 303   // Fix references to daughter tracks\r
 304   //  if(fFixRefs) {\r
 305   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();\r
 306   //    fixer->FixReferences((AliAODEvent*)event);\r
 307   //    delete fixer;\r
 308   //  }\r
 309   //\r
 310 \r
 311 \r
 312   // physics selection requirements\r
 313   if(fUsePhysicsSelection){\r
 314     Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);\r
 315     if(!isSelected) {\r
 316       if(accept) fWhyRejection=7;\r
 317       fEvRejectionBits+=1<<kPhysicsSelection;\r
 318       accept=kFALSE;\r
 319     }else{\r
 320       if(fUseOnlyOneTrigger){\r
 321         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){\r
 322           if(accept) fWhyRejection=7;\r
 323           fEvRejectionBits+=1<<kPhysicsSelection;\r
 324           accept=kFALSE;\r
 325         }\r
 326       }\r
 327     }\r
 328   }\r
 329 \r
 330   // vertex requirements\r
 331    \r
 332   const AliVVertex *vertex = event->GetPrimaryVertex();\r
 333 \r
 334   if(!vertex){\r
 335     accept=kFALSE;\r
 336     fEvRejectionBits+=1<<kNoVertex;\r
 337   }else{\r
 338     TString title=vertex->GetTitle();\r
 339     if(title.Contains("Z") && fMinVtxType>1){\r
 340       accept=kFALSE;\r
 341       fEvRejectionBits+=1<<kNoVertex;\r
 342     }\r
 343     else if(title.Contains("3D") && fMinVtxType>2){\r
 344       accept=kFALSE;\r
 345       fEvRejectionBits+=1<<kNoVertex;\r
 346     }\r
 347     if(vertex->GetNContributors()<fMinVtxContr){\r
 348       accept=kFALSE;\r
 349       fEvRejectionBits+=1<<kTooFewVtxContrib;\r
 350     }\r
 351     if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {\r
 352       fEvRejectionBits+=1<<kZVtxOutFid;\r
 353       if(accept) fWhyRejection=6;\r
 354       accept=kFALSE;\r
 355     } \r
 356   }\r
 357 \r
 358 \r
 359   // pile-up rejection\r
 360   if(fOptPileup==kRejectPileupEvent){\r
 361     Int_t cutc=(Int_t)fMinContrPileup;\r
 362     Double_t cutz=(Double_t)fMinDzPileup;\r
 363     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {\r
 364       if(accept) fWhyRejection=1;\r
 365       fEvRejectionBits+=1<<kPileupSPD;\r
 366       accept=kFALSE;\r
 367     }\r
 368   }\r
 369 \r
 370   // centrality selection\r
 371   if (fUseCentrality!=kCentOff) {  \r
 372     Int_t rejection=IsEventSelectedInCentrality(event);    \r
 373     if(rejection>1){      \r
 374       if(accept) fWhyRejection=rejection;      \r
 375       fEvRejectionBits+=1<<kOutsideCentrality;\r
 376       accept=kFALSE;\r
 377     }\r
 378   }\r
 379 \r
 380   return accept;\r
 381 }\r
 382 //---------------------------------------------------------------------------\r
 383 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {\r
 384   //\r
 385   // Daughter tracks selection\r
 386   // \r
 387   if(!fTrackCuts) return kTRUE;\r
 388  \r
 389   Int_t ndaughters = d->GetNDaughters();\r
 390   AliAODVertex *vAOD = d->GetPrimaryVtx();\r
 391   Double_t pos[3],cov[6];\r
 392   vAOD->GetXYZ(pos);\r
 393   vAOD->GetCovarianceMatrix(cov);\r
 394   const AliESDVertex vESD(pos,cov,100.,100);\r
 395   \r
 396   Bool_t retval=kTRUE;\r
 397   \r
 398   for(Int_t idg=0; idg<ndaughters; idg++) {\r
 399     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);\r
 400     if(!dgTrack) {retval = kFALSE; continue;}\r
 401     //printf("charge %d\n",dgTrack->Charge());\r
 402     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0\r
 403 \r
 404     if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)\r
 405       { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }\r
 406 \r
 407     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;\r
 408   }\r
 409   \r
 410   return retval;\r
 411 }\r
 412 //---------------------------------------------------------------------------\r
 413 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {\r
 414   //\r
 415   // Convert to ESDtrack, relate to vertex and check cuts\r
 416   //\r
 417   if(!cuts) return kTRUE;\r
 418 \r
 419   Bool_t retval=kTRUE;\r
 420 \r
 421   // convert to ESD track here\r
 422   AliESDtrack esdTrack(track);\r
 423   // set the TPC cluster info\r
 424   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());\r
 425   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());\r
 426   esdTrack.SetTPCPointsF(track->GetTPCNclsF());\r
 427   // needed to calculate the impact parameters\r
 428   esdTrack.RelateToVertex(primary,0.,3.); \r
 429 \r
 430   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;\r
 431  \r
 432   if(fOptPileup==kRejectTracksFromPileupVertex){\r
 433     // to be implemented\r
 434     // we need either to have here the AOD Event, \r
 435     // or to have the pileup vertex object\r
 436   }\r
 437   return retval; \r
 438 }\r
 439 //---------------------------------------------------------------------------\r
 440 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {\r
 441   // Set the pt bins\r
 442 \r
 443   if(fPtBinLimits) {\r
 444     delete [] fPtBinLimits;\r
 445     fPtBinLimits = NULL;\r
 446     printf("Changing the pt bins\n");\r
 447   }\r
 448 \r
 449   if(nPtBinLimits != fnPtBins+1){\r
 450     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;\r
 451     SetNPtBins(nPtBinLimits-1);\r
 452   }\r
 453 \r
 454   fnPtBinLimits = nPtBinLimits;\r
 455   SetGlobalIndex();\r
 456   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;\r
 457   fPtBinLimits = new Float_t[fnPtBinLimits];\r
 458   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];\r
 459 \r
 460   return;\r
 461 }\r
 462 //---------------------------------------------------------------------------\r
 463 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){\r
 464   // Set the variable names\r
 465 \r
 466   if(fVarNames) {\r
 467     delete [] fVarNames;\r
 468     fVarNames = NULL;\r
 469     //printf("Changing the variable names\n");\r
 470   }\r
 471   if(nVars!=fnVars){\r
 472     printf("Wrong number of variables: it has to be %d\n",fnVars);\r
 473     return;\r
 474   }\r
 475   //fnVars=nVars;\r
 476   fVarNames = new TString[nVars];\r
 477   fIsUpperCut = new Bool_t[nVars];\r
 478   for(Int_t iv=0; iv<nVars; iv++) {\r
 479     fVarNames[iv] = varNames[iv];\r
 480     fIsUpperCut[iv] = isUpperCut[iv];\r
 481   }\r
 482 \r
 483   return;\r
 484 }\r
 485 //---------------------------------------------------------------------------\r
 486 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {\r
 487   // Set the variables to be used for cuts optimization\r
 488 \r
 489   if(fVarsForOpt) {\r
 490     delete [] fVarsForOpt;\r
 491     fVarsForOpt = NULL;\r
 492     //printf("Changing the variables for cut optimization\n");\r
 493   }\r
 494   \r
 495   if(nVars==0){//!=fnVars) {\r
 496     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);\r
 497     return;\r
 498   } \r
 499   \r
 500   fnVarsForOpt = 0;\r
 501   fVarsForOpt = new Bool_t[fnVars];\r
 502   for(Int_t iv=0; iv<fnVars; iv++) {\r
 503     fVarsForOpt[iv]=forOpt[iv];\r
 504     if(fVarsForOpt[iv]) fnVarsForOpt++;\r
 505   }\r
 506 \r
 507   return;\r
 508 }\r
 509 \r
 510 //---------------------------------------------------------------------------\r
 511 void AliRDHFCuts::SetUseCentrality(Int_t flag) {\r
 512   //\r
 513   // set centrality estimator  \r
 514   //\r
 515   fUseCentrality=flag;\r
 516   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");\r
 517  \r
 518   return;\r
 519 }\r
 520 \r
 521 \r
 522 //---------------------------------------------------------------------------\r
 523 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {\r
 524   //\r
 525   // store the cuts\r
 526   //\r
 527   if(nVars!=fnVars) {\r
 528     printf("Wrong number of variables: it has to be %d\n",fnVars);\r
 529     AliFatal("exiting");\r
 530   } \r
 531   if(nPtBins!=fnPtBins) {\r
 532     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);\r
 533     AliFatal("exiting");\r
 534   } \r
 535 \r
 536   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];\r
 537   \r
 538 \r
 539   for(Int_t iv=0; iv<fnVars; iv++) {\r
 540 \r
 541     for(Int_t ib=0; ib<fnPtBins; ib++) {\r
 542 \r
 543       //check\r
 544       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {\r
 545         cout<<"Overflow, exit..."<<endl;\r
 546         return;\r
 547       }\r
 548 \r
 549       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];\r
 550 \r
 551     }\r
 552   }\r
 553   return;\r
 554 }\r
 555 //---------------------------------------------------------------------------\r
 556 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){\r
 557   //\r
 558   // store the cuts\r
 559   //\r
 560   if(glIndex != fGlobalIndex){\r
 561     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;\r
 562     AliFatal("exiting");\r
 563   }\r
 564   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];\r
 565 \r
 566   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){\r
 567     fCutsRD[iGl] = cutsRDGlob[iGl];\r
 568   }\r
 569   return;\r
 570 }\r
 571 //---------------------------------------------------------------------------\r
 572 void AliRDHFCuts::PrintAll() const {\r
 573   //\r
 574   // print all cuts values\r
 575   // \r
 576 \r
 577   printf("Minimum vtx type %d\n",fMinVtxType);\r
 578   printf("Minimum vtx contr %d\n",fMinVtxContr);\r
 579   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);\r
 580   printf("Min SPD mult %d\n",fMinSPDMultiplicity);\r
 581   printf("Use PID %d\n",(Int_t)fUsePID);\r
 582   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);\r
 583   printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");\r
 584   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");\r
 585   if(fOptPileup==1) printf(" -- Reject pileup event");\r
 586   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");\r
 587   if(fUseCentrality>0) {\r
 588     TString estimator="";\r
 589     if(fUseCentrality==1) estimator = "V0";\r
 590     if(fUseCentrality==2) estimator = "Tracks";\r
 591     if(fUseCentrality==3) estimator = "Tracklets";\r
 592     if(fUseCentrality==4) estimator = "SPD clusters outer"; \r
 593     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());\r
 594   }\r
 595   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);\r
 596 \r
 597   if(fVarNames){\r
 598     cout<<"Array of variables"<<endl;\r
 599     for(Int_t iv=0;iv<fnVars;iv++){\r
 600       cout<<fVarNames[iv]<<"\t";\r
 601     }\r
 602     cout<<endl;\r
 603   }\r
 604   if(fVarsForOpt){\r
 605     cout<<"Array of optimization"<<endl;\r
 606     for(Int_t iv=0;iv<fnVars;iv++){\r
 607       cout<<fVarsForOpt[iv]<<"\t";\r
 608     }\r
 609     cout<<endl;\r
 610   }\r
 611   if(fIsUpperCut){\r
 612     cout<<"Array of upper/lower cut"<<endl;\r
 613    for(Int_t iv=0;iv<fnVars;iv++){\r
 614      cout<<fIsUpperCut[iv]<<"\t";\r
 615    }\r
 616    cout<<endl;\r
 617   }\r
 618   if(fPtBinLimits){\r
 619     cout<<"Array of ptbin limits"<<endl;\r
 620     for(Int_t ib=0;ib<fnPtBinLimits;ib++){\r
 621       cout<<fPtBinLimits[ib]<<"\t";\r
 622     }\r
 623     cout<<endl;\r
 624   }\r
 625   if(fCutsRD){\r
 626     cout<<"Matrix of cuts"<<endl;\r
 627    for(Int_t iv=0;iv<fnVars;iv++){\r
 628      for(Int_t ib=0;ib<fnPtBins;ib++){\r
 629        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";\r
 630      } \r
 631      cout<<endl;\r
 632    }\r
 633    cout<<endl;\r
 634   }\r
 635   return;\r
 636 }\r
 637 //---------------------------------------------------------------------------\r
 638 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{\r
 639   //\r
 640   // get the cuts\r
 641   //\r
 642 \r
 643   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;\r
 644 \r
 645 \r
 646   Int_t iv,ib;\r
 647   if(!cutsRD) {\r
 648     //cout<<"Initialization..."<<endl;\r
 649     cutsRD=new Float_t*[fnVars];\r
 650     for(iv=0; iv<fnVars; iv++) {\r
 651       cutsRD[iv] = new Float_t[fnPtBins];\r
 652     }\r
 653   }\r
 654   \r
 655   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {\r
 656     GetVarPtIndex(iGlobal,iv,ib);\r
 657     cutsRD[iv][ib] = fCutsRD[iGlobal];\r
 658   }\r
 659 \r
 660   return;\r
 661 }\r
 662 \r
 663 //---------------------------------------------------------------------------\r
 664 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{\r
 665   //\r
 666   // give the global index from variable and pt bin\r
 667   //\r
 668   return iPtBin*fnVars+iVar;\r
 669 }\r
 670 \r
 671 //---------------------------------------------------------------------------\r
 672 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {\r
 673   //\r
 674   //give the index of the variable and of the pt bin from the global index\r
 675   //\r
 676   iPtBin=(Int_t)iGlob/fnVars;\r
 677   iVar=iGlob%fnVars;\r
 678 \r
 679   return;\r
 680 }\r
 681 \r
 682 //---------------------------------------------------------------------------\r
 683 Int_t AliRDHFCuts::PtBin(Double_t pt) const {\r
 684   //\r
 685   //give the pt bin where the pt lies.\r
 686   //\r
 687   Int_t ptbin=-1;\r
 688   if(pt<fPtBinLimits[0])return ptbin;\r
 689   for (Int_t i=0;i<fnPtBins;i++){\r
 690     if(pt<fPtBinLimits[i+1]) {\r
 691       ptbin=i;\r
 692       break;\r
 693     }\r
 694   }\r
 695   return ptbin;\r
 696 }\r
 697 //-------------------------------------------------------------------\r
 698 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {\r
 699   // \r
 700   // Give the value of cut set for the variable iVar and the pt bin iPtBin\r
 701   //\r
 702   if(!fCutsRD){\r
 703     cout<<"Cuts not iniziaisez yet"<<endl;\r
 704     return 0;\r
 705   }\r
 706   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];\r
 707 }\r
 708 //-------------------------------------------------------------------\r
 709 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {\r
 710   //\r
 711   // Get centrality percentile\r
 712   //\r
 713 \r
 714   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 715   if(mcArray) {fUseAOD049=kFALSE;}\r
 716 \r
 717   AliAODHeader *header=aodEvent->GetHeader();\r
 718   AliCentrality *centrality=header->GetCentralityP();\r
 719   Float_t cent=-999.;\r
 720   Bool_t isSelRun=kFALSE;\r
 721   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};\r
 722   if(!centrality) return cent;\r
 723   else{\r
 724     if (estimator==kCentV0M){\r
 725       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));\r
 726       if(cent<0){\r
 727         Int_t quality = centrality->GetQuality();\r
 728         if(quality<=1){\r
 729           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
 730         }else{\r
 731           Int_t runnum=aodEvent->GetRunNumber();\r
 732           for(Int_t ir=0;ir<5;ir++){\r
 733             if(runnum==selRun[ir]){\r
 734               isSelRun=kTRUE;\r
 735               break;\r
 736             }\r
 737           }\r
 738           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
 739         }\r
 740       }\r
 741 \r
 742       //temporary fix for AOD049 outliers\r
 743       if(fUseAOD049&&cent>=0){\r
 744         Float_t v0=0;\r
 745         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();\r
 746         v0+=aodV0->GetMTotV0A();\r
 747         v0+=aodV0->GetMTotV0C();\r
 748         if(cent==0&&v0<19500)return -1;//filtering issue\r
 749         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());\r
 750         Float_t val= 1.30552 +  0.147931 * v0;\r
 751         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};\r
 752         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier\r
 753       }\r
 754     }\r
 755     else {\r
 756       if (estimator==kCentTRK) {\r
 757         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));\r
 758         if(cent<0){\r
 759           Int_t quality = centrality->GetQuality();\r
 760           if(quality<=1){\r
 761             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
 762           }else{\r
 763             Int_t runnum=aodEvent->GetRunNumber();\r
 764             for(Int_t ir=0;ir<5;ir++){\r
 765               if(runnum==selRun[ir]){\r
 766                 isSelRun=kTRUE;\r
 767                 break;\r
 768               }\r
 769             }\r
 770             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
 771           }\r
 772         }\r
 773       }\r
 774       else{\r
 775         if (estimator==kCentTKL){\r
 776           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));\r
 777           if(cent<0){\r
 778             Int_t quality = centrality->GetQuality();\r
 779             if(quality<=1){\r
 780               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
 781             }else{\r
 782               Int_t runnum=aodEvent->GetRunNumber();\r
 783               for(Int_t ir=0;ir<5;ir++){\r
 784                 if(runnum==selRun[ir]){\r
 785                   isSelRun=kTRUE;\r
 786                   break;\r
 787             }\r
 788               }\r
 789               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
 790             }   \r
 791           }\r
 792         }\r
 793         else{\r
 794           if (estimator==kCentCL1){\r
 795             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));\r
 796             if(cent<0){\r
 797               Int_t quality = centrality->GetQuality();\r
 798               if(quality<=1){\r
 799                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
 800               }else{\r
 801                 Int_t runnum=aodEvent->GetRunNumber();\r
 802                 for(Int_t ir=0;ir<5;ir++){\r
 803                   if(runnum==selRun[ir]){\r
 804                     isSelRun=kTRUE;\r
 805                     break;\r
 806                   }\r
 807                 }\r
 808                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
 809               }\r
 810             }\r
 811           }\r
 812           else {\r
 813             AliWarning("Centrality estimator not valid");\r
 814             \r
 815           }\r
 816         }\r
 817       }\r
 818     } \r
 819   }\r
 820   return cent;\r
 821 }\r
 822 //-------------------------------------------------------------------\r
 823 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {\r
 824   //\r
 825   // Compare two cuts objects\r
 826   //\r
 827 \r
 828   Bool_t areEqual=kTRUE;\r
 829 \r
 830   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}\r
 831 \r
 832   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}\r
 833 \r
 834   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}\r
 835 \r
 836   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}\r
 837 \r
 838   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}\r
 839 \r
 840   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}\r
 841   if(fTrackCuts){\r
 842     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}\r
 843 \r
 844     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}\r
 845 \r
 846     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}\r
 847 \r
 848     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;}\r
 849   }\r
 850 \r
 851   if(fCutsRD) {\r
 852    for(Int_t iv=0;iv<fnVars;iv++) {\r
 853      for(Int_t ib=0;ib<fnPtBins;ib++) {\r
 854        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {\r
 855          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";\r
 856          areEqual=kFALSE;\r
 857        }\r
 858      }\r
 859    }\r
 860   }\r
 861 \r
 862   return areEqual;\r
 863 }\r
 864 //---------------------------------------------------------------------------\r
 865 void AliRDHFCuts::MakeTable() const {\r
 866   //\r
 867   // print cuts values in table format\r
 868   // \r
 869 \r
 870         TString ptString = "pT range";\r
 871         if(fVarNames && fPtBinLimits && fCutsRD){\r
 872                 TString firstLine(Form("*       %-15s",ptString.Data()));\r
 873                 for (Int_t ivar=0; ivar<fnVars; ivar++){\r
 874                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());\r
 875                         if (ivar == fnVars){\r
 876                                 firstLine+="*\n";\r
 877                         }\r
 878                 }\r
 879                 Printf("%s",firstLine.Data());\r
 880                 \r
 881                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){\r
 882                         TString line;\r
 883                         if (ipt==fnPtBins-1){\r
 884                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);\r
 885                         }\r
 886                         else{\r
 887                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);\r
 888                         }\r
 889                         for (Int_t ivar=0; ivar<fnVars; ivar++){\r
 890                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);\r
 891                         }\r
 892                         Printf("%s",line.Data());\r
 893                 }\r
 894 \r
 895         }\r
 896 \r
 897 \r
 898   return;\r
 899 }\r
 900 //--------------------------------------------------------------------------\r
 901 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
 902                                         AliAODEvent *aod) const\r
 903 {\r
 904   //\r
 905   // Recalculate primary vertex without daughters\r
 906   //\r
 907 \r
 908   if(!aod) {\r
 909     AliError("Can not remove daughters from vertex without AOD event");\r
 910     return 0;\r
 911   }   \r
 912 \r
 913   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);\r
 914   if(!recvtx){\r
 915     AliDebug(2,"Removal of daughter tracks failed");\r
 916     return kFALSE;\r
 917   }\r
 918 \r
 919 \r
 920   //set recalculed primary vertex\r
 921   d->SetOwnPrimaryVtx(recvtx);\r
 922   delete recvtx;\r
 923 \r
 924   return kTRUE;\r
 925 }\r
 926 //--------------------------------------------------------------------------\r
 927 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const\r
 928 {\r
 929   //\r
 930   // Recalculate primary vertex without daughters\r
 931   //\r
 932 \r
 933   if(!aod) {\r
 934     AliError("Can not get MC vertex without AOD event");\r
 935     return kFALSE;\r
 936   }   \r
 937 \r
 938   // load MC header\r
 939   AliAODMCHeader *mcHeader = \r
 940     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
 941   if(!mcHeader) {\r
 942     AliError("Can not get MC vertex without AODMCHeader event");\r
 943     return kFALSE;\r
 944   }\r
 945   Double_t pos[3];\r
 946   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};\r
 947   mcHeader->GetVertex(pos);\r
 948   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);\r
 949 \r
 950   if(!recvtx){\r
 951     AliDebug(2,"Removal of daughter tracks failed");\r
 952     return kFALSE;\r
 953   }\r
 954 \r
 955   //set recalculed primary vertex\r
 956   d->SetOwnPrimaryVtx(recvtx);\r
 957 \r
 958   d->RecalculateImpPars(recvtx,aod);\r
 959 \r
 960   delete recvtx;\r
 961 \r
 962   return kTRUE;\r
 963 }\r
 964 //--------------------------------------------------------------------------\r
 965 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
 966                                      AliAODEvent *aod,\r
 967                                      AliAODVertex *origownvtx) const\r
 968 {\r
 969   //\r
 970   // Clean-up own primary vertex if needed\r
 971   //\r
 972 \r
 973   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {\r
 974     d->UnsetOwnPrimaryVtx();\r
 975     if(origownvtx) {\r
 976       d->SetOwnPrimaryVtx(origownvtx);\r
 977       delete origownvtx; origownvtx=NULL;\r
 978     }\r
 979     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);\r
 980   } else {\r
 981     if(origownvtx) {\r
 982       delete origownvtx; origownvtx=NULL;\r
 983     }\r
 984   }\r
 985   return;\r
 986 }\r
 987 //--------------------------------------------------------------------------\r
 988 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const \r
 989 {\r
 990   //\r
 991   // Checks if this candidate is matched to MC signal\r
 992   // \r
 993 \r
 994   if(!aod) return kFALSE;\r
 995 \r
 996   // get the MC array\r
 997   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 998 \r
 999   if(!mcArray) return kFALSE;\r
1000 \r
1001   // try to match  \r
1002   Int_t label = d->MatchToMC(pdg,mcArray);\r
1003   \r
1004   if(label>=0) {\r
1005     //printf("MATCH!\n");\r
1006     return kTRUE;\r
1007   }\r
1008 \r
1009   return kFALSE;\r
1010 }\r
1011 \r
1012 \r