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