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
  47 using std::cout;
  48 using std::endl;
  49
  50 ClassImp(AliRDHFCuts)
  51
  52 //--------------------------------------------------------------------------
  53 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
  54 AliAnalysisCuts(name,title),
  55 fMinVtxType(3),
  56 fMinVtxContr(1),
  57 fMaxVtxRedChi2(1e6),
  58 fMaxVtxZ(10.),
  59 fMinSPDMultiplicity(0),
  60 fTriggerMask(AliVEvent::kAnyINT),
  61 fUseOnlyOneTrigger(kFALSE),
  62 fTrackCuts(0),
  63 fnPtBins(1),
  64 fnPtBinLimits(1),
  65 fPtBinLimits(0),
  66 fnVars(1),
  67 fVarNames(0),
  68 fnVarsForOpt(0),
  69 fVarsForOpt(0),
  70 fGlobalIndex(1),
  71 fCutsRD(0),
  72 fIsUpperCut(0),
  73 fUsePID(kFALSE),
  74 fUseAOD049(kFALSE),
  75 fPidHF(0),
  76 fWhyRejection(0),
  77 fEvRejectionBits(0),
  78 fRemoveDaughtersFromPrimary(kFALSE),
  79 fRecomputePrimVertex(kFALSE),
  80 fUseMCVertex(kFALSE),
  81 fUsePhysicsSelection(kTRUE),
  82 fOptPileup(0),
  83 fMinContrPileup(3),
  84 fMinDzPileup(0.6),
  85 fUseCentrality(0),
  86 fMinCentrality(0.),
  87 fMaxCentrality(100.),
  88 fFixRefs(kFALSE),
  89 fIsSelectedCuts(0),
  90 fIsSelectedPID(0),
  91 fMinPtCand(-1.),
  92 fMaxPtCand(100000.),
  93 fKeepSignalMC(kFALSE),
  94 fIsCandTrackSPDFirst(kFALSE),
  95 fMaxPtCandTrackSPDFirst(0.),
  96 fApplySPDDeadPbPb2011(kFALSE),
  97 fRemoveTrackletOutliers(kFALSE),
  98 fKinkReject(kFALSE)
  99 {
 100   //
 101   // Default Constructor
 102   //
 103   fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
 104 }
 105 //--------------------------------------------------------------------------
 106 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
 107   AliAnalysisCuts(source),
 108   fMinVtxType(source.fMinVtxType),
 109   fMinVtxContr(source.fMinVtxContr),
 110   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
 111   fMaxVtxZ(source.fMaxVtxZ),
 112   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
 113   fTriggerMask(source.fTriggerMask),
 114   fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
 115   fTriggerClass(),
 116   fTrackCuts(0),
 117   fnPtBins(source.fnPtBins),
 118   fnPtBinLimits(source.fnPtBinLimits),
 119   fPtBinLimits(0),
 120   fnVars(source.fnVars),
 121   fVarNames(0),
 122   fnVarsForOpt(source.fnVarsForOpt),
 123   fVarsForOpt(0),
 124   fGlobalIndex(source.fGlobalIndex),
 125   fCutsRD(0),
 126   fIsUpperCut(0),
 127   fUsePID(source.fUsePID),
 128   fUseAOD049(source.fUseAOD049),
 129   fPidHF(0),
 130   fWhyRejection(source.fWhyRejection),
 131   fEvRejectionBits(source.fEvRejectionBits),
 132   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
 133   fRecomputePrimVertex(source.fRecomputePrimVertex),
 134   fUseMCVertex(source.fUseMCVertex),
 135   fUsePhysicsSelection(source.fUsePhysicsSelection),
 136   fOptPileup(source.fOptPileup),
 137   fMinContrPileup(source.fMinContrPileup),
 138   fMinDzPileup(source.fMinDzPileup),
 139   fUseCentrality(source.fUseCentrality),
 140   fMinCentrality(source.fMinCentrality),
 141   fMaxCentrality(source.fMaxCentrality),
 142   fFixRefs(source.fFixRefs),
 143   fIsSelectedCuts(source.fIsSelectedCuts),
 144   fIsSelectedPID(source.fIsSelectedPID),
 145   fMinPtCand(source.fMinPtCand),
 146   fMaxPtCand(source.fMaxPtCand),
 147   fKeepSignalMC(source.fKeepSignalMC),
 148   fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
 149   fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
 150   fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
 151   fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
 152   fKinkReject(source.fKinkReject)
 153 {
 154   //
 155   // Copy constructor
 156   //
 157   cout<<"Copy constructor"<<endl;
 158   fTriggerClass[0] = source.fTriggerClass[0];
 159   fTriggerClass[1] = source.fTriggerClass[1];
 160   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
 161   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 162   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 163   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 164   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 165   if(source.fPidHF) SetPidHF(source.fPidHF);
 166   PrintAll();
 167
 168 }
 169 //--------------------------------------------------------------------------
 170 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
 171 {
 172   //
 173   // assignment operator
 174   //
 175   if(&source == this) return *this;
 176
 177   AliAnalysisCuts::operator=(source);
 178
 179   fMinVtxType=source.fMinVtxType;
 180   fMinVtxContr=source.fMinVtxContr;
 181   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
 182   fMaxVtxZ=source.fMaxVtxZ;
 183   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
 184   fTriggerMask=source.fTriggerMask;
 185   fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
 186   fTriggerClass[0]=source.fTriggerClass[0];
 187   fTriggerClass[1]=source.fTriggerClass[1];
 188   fnPtBins=source.fnPtBins;
 189   fnPtBinLimits=source.fnPtBinLimits;
 190   fnVars=source.fnVars;
 191   fGlobalIndex=source.fGlobalIndex;
 192   fnVarsForOpt=source.fnVarsForOpt;
 193   fUsePID=source.fUsePID;
 194   fUseAOD049=source.fUseAOD049;
 195   if(fPidHF) delete fPidHF;
 196   fPidHF=new AliAODPidHF(*(source.GetPidHF()));
 197   fWhyRejection=source.fWhyRejection;
 198   fEvRejectionBits=source.fEvRejectionBits;
 199   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
 200   fRecomputePrimVertex=source.fRecomputePrimVertex;
 201   fUseMCVertex=source.fUseMCVertex;
 202   fUsePhysicsSelection=source.fUsePhysicsSelection;
 203   fOptPileup=source.fOptPileup;
 204   fMinContrPileup=source.fMinContrPileup;
 205   fMinDzPileup=source.fMinDzPileup;
 206   fUseCentrality=source.fUseCentrality;
 207   fMinCentrality=source.fMinCentrality;
 208   fMaxCentrality=source.fMaxCentrality;
 209   fFixRefs=source.fFixRefs;
 210   fIsSelectedCuts=source.fIsSelectedCuts;
 211   fIsSelectedPID=source.fIsSelectedPID;
 212   fMinPtCand=source.fMinPtCand;
 213   fMaxPtCand=source.fMaxPtCand;
 214   fKeepSignalMC=source.fKeepSignalMC;
 215   fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
 216   fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
 217   fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
 218   fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
 219   fKinkReject=source.fKinkReject;
 220
 221   if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
 222   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 223   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 224   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 225   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 226   PrintAll();
 227
 228   return *this;
 229 }
 230 //--------------------------------------------------------------------------
 231 AliRDHFCuts::~AliRDHFCuts() {
 232   //
 233   // Default Destructor
 234   //
 235   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
 236   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
 237   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
 238   if(fCutsRD) {
 239     delete [] fCutsRD;
 240     fCutsRD=0;
 241   }
 242   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
 243   if(fPidHF){
 244     delete fPidHF;
 245     fPidHF=0;
 246   }
 247 }
 248 //---------------------------------------------------------------------------
 249 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
 250   //
 251   // Centrality selection
 252   //
 253   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
 254     AliWarning("Centrality estimator not valid");
 255     return 3;
 256   }else{
 257     Float_t centvalue=GetCentrality((AliAODEvent*)event);
 258     if (centvalue<-998.){//-999 if no centralityP
 259       return 0;
 260     }else{
 261       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
 262         return 2;
 263       }
 264     }
 265   }
 266   return 0;
 267 }
 268 //---------------------------------------------------------------------------
 269 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
 270   //
 271   // Event selection
 272   //
 273   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
 274
 275   fWhyRejection=0;
 276   fEvRejectionBits=0;
 277   Bool_t accept=kTRUE;
 278
 279   if(fRecomputePrimVertex){
 280     Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
 281     if(!vertOK){
 282       fWhyRejection=6;
 283       return kFALSE;
 284     }
 285   }
 286
 287   // check if it's MC
 288   Bool_t isMC=kFALSE;
 289   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 290   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
 291
 292   // settings for the TPC dE/dx BB parameterization
 293   if(fPidHF){
 294     if(fPidHF->GetPidResponse()==0x0){
 295       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 296       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 297       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 298       fPidHF->SetPidResponse(pidResp);
 299     }
 300     if(fPidHF->GetOldPid()) {
 301       // pp, from LHC10d onwards
 302       if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
 303          event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
 304       // pp, 2011 low energy run
 305       if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
 306         fPidHF->SetppLowEn2011(kTRUE);
 307         fPidHF->SetOnePad(kFALSE);
 308       }
 309       // PbPb LHC10h
 310       if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
 311       // MC
 312       if(isMC) fPidHF->SetMC(kTRUE);
 313       if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
 314         fPidHF->SetMClowenpp2011(kTRUE);
 315       fPidHF->SetBetheBloch();
 316     }else{
 317       // check that AliPIDResponse object was properly set in case of using OADB
 318       if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
 319     }
 320   }
 321
 322
 323   // trigger class
 324   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
 325   // don't do for MC and for PbPb 2010 data
 326   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
 327     if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
 328        (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
 329       fWhyRejection=5;
 330       fEvRejectionBits+=1<<kNotSelTrigger;
 331       accept=kFALSE;
 332     }
 333   }
 334
 335   // TEMPORARY FIX FOR REFERENCES
 336   // Fix references to daughter tracks
 337   //  if(fFixRefs) {
 338   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
 339   //    fixer->FixReferences((AliAODEvent*)event);
 340   //    delete fixer;
 341   //  }
 342   //
 343
 344
 345   // physics selection requirements
 346   if(fUsePhysicsSelection){
 347     Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
 348     if(!isSelected) {
 349       if(accept) fWhyRejection=7;
 350       fEvRejectionBits+=1<<kPhysicsSelection;
 351       accept=kFALSE;
 352     }else{
 353       if(fUseOnlyOneTrigger){
 354         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
 355           if(accept) fWhyRejection=7;
 356           fEvRejectionBits+=1<<kPhysicsSelection;
 357           accept=kFALSE;
 358         }
 359       }
 360     }
 361   }
 362
 363   // vertex requirements
 364
 365   const AliVVertex *vertex = event->GetPrimaryVertex();
 366
 367   if(!vertex){
 368     accept=kFALSE;
 369     fEvRejectionBits+=1<<kNoVertex;
 370   }else{
 371     TString title=vertex->GetTitle();
 372     if(title.Contains("Z") && fMinVtxType>1){
 373       accept=kFALSE;
 374       fEvRejectionBits+=1<<kNoVertex;
 375     }
 376     else if(title.Contains("3D") && fMinVtxType>2){
 377       accept=kFALSE;
 378       fEvRejectionBits+=1<<kNoVertex;
 379     }
 380     if(vertex->GetNContributors()<fMinVtxContr){
 381       accept=kFALSE;
 382       fEvRejectionBits+=1<<kTooFewVtxContrib;
 383     }
 384     if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
 385       fEvRejectionBits+=1<<kZVtxOutFid;
 386       if(accept) fWhyRejection=6;
 387       accept=kFALSE;
 388     }
 389   }
 390
 391
 392   // pile-up rejection
 393   if(fOptPileup==kRejectPileupEvent){
 394     Int_t cutc=(Int_t)fMinContrPileup;
 395     Double_t cutz=(Double_t)fMinDzPileup;
 396     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
 397       if(accept) fWhyRejection=1;
 398       fEvRejectionBits+=1<<kPileupSPD;
 399       accept=kFALSE;
 400     }
 401   }
 402
 403   // centrality selection
 404   if (fUseCentrality!=kCentOff) {
 405     Int_t rejection=IsEventSelectedInCentrality(event);
 406     if(rejection>1){
 407       if(accept) fWhyRejection=rejection;
 408       fEvRejectionBits+=1<<kOutsideCentrality;
 409       accept=kFALSE;
 410     }
 411   }
 412
 413   // PbPb2011 outliers in tracklets vs. VZERO
 414   if(!isMC && event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
 415     if(fRemoveTrackletOutliers){
 416       Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
 417       Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
 418       Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
 419       if(v0cent<cutval){
 420         if(accept) fWhyRejection=2;
 421         fEvRejectionBits+=1<<kOutsideCentrality;
 422          accept=kFALSE;
 423       }
 424     }
 425   }
 426
 427   return accept;
 428 }
 429 //---------------------------------------------------------------------------
 430 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
 431   //
 432   // Daughter tracks selection
 433   //
 434   if(!fTrackCuts) return kTRUE;
 435
 436   Int_t ndaughters = d->GetNDaughters();
 437   AliAODVertex *vAOD = d->GetPrimaryVtx();
 438   Double_t pos[3],cov[6];
 439   vAOD->GetXYZ(pos);
 440   vAOD->GetCovarianceMatrix(cov);
 441   const AliESDVertex vESD(pos,cov,100.,100);
 442
 443   Bool_t retval=kTRUE;
 444
 445   for(Int_t idg=0; idg<ndaughters; idg++) {
 446     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
 447     if(!dgTrack) {retval = kFALSE; continue;}
 448     //printf("charge %d\n",dgTrack->Charge());
 449     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
 450
 451     if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
 452       { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
 453
 454     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
 455   }
 456
 457   return retval;
 458 }
 459 //---------------------------------------------------------------------------
 460 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
 461   //
 462   // Convert to ESDtrack, relate to vertex and check cuts
 463   //
 464   if(!cuts) return kTRUE;
 465
 466   if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
 467
 468   Bool_t retval=kTRUE;
 469
 470   // convert to ESD track here
 471   AliESDtrack esdTrack(track);
 472   // set the TPC cluster info
 473   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
 474   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
 475   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
 476   // needed to calculate the impact parameters
 477   esdTrack.RelateToVertex(primary,0.,3.);
 478
 479   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
 480
 481
 482   if(fKinkReject){
 483    AliAODVertex *maybeKink=track->GetProdVertex();
 484    if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
 485   }
 486
 487   if(fOptPileup==kRejectTracksFromPileupVertex){
 488     // to be implemented
 489     // we need either to have here the AOD Event,
 490     // or to have the pileup vertex object
 491   }
 492
 493   if(fApplySPDDeadPbPb2011){
 494     Bool_t deadSPDLay1PbPb2011[20][4]={
 495       {kTRUE,kTRUE,kTRUE,kTRUE},
 496       {kTRUE,kTRUE,kTRUE,kTRUE},
 497       {kTRUE,kTRUE,kTRUE,kTRUE},
 498       {kTRUE,kTRUE,kTRUE,kTRUE},
 499       {kTRUE,kTRUE,kTRUE,kTRUE},
 500       {kFALSE,kFALSE,kTRUE,kTRUE},
 501       {kTRUE,kTRUE,kFALSE,kFALSE},
 502       {kTRUE,kTRUE,kTRUE,kTRUE},
 503       {kFALSE,kFALSE,kTRUE,kTRUE},
 504       {kTRUE,kTRUE,kTRUE,kTRUE},
 505       {kTRUE,kTRUE,kFALSE,kFALSE},
 506       {kTRUE,kTRUE,kTRUE,kTRUE},
 507       {kFALSE,kFALSE,kFALSE,kFALSE},
 508       {kFALSE,kFALSE,kTRUE,kTRUE},
 509       {kFALSE,kFALSE,kFALSE,kFALSE},
 510       {kFALSE,kFALSE,kFALSE,kFALSE},
 511       {kTRUE,kTRUE,kTRUE,kTRUE},
 512       {kTRUE,kTRUE,kFALSE,kFALSE},
 513       {kFALSE,kFALSE,kFALSE,kFALSE},
 514       {kFALSE,kFALSE,kFALSE,kFALSE}
 515     };
 516     Bool_t deadSPDLay2PbPb2011[40][4]={
 517       {kTRUE,kTRUE,kTRUE,kTRUE},
 518       {kTRUE,kTRUE,kTRUE,kTRUE},
 519       {kTRUE,kTRUE,kTRUE,kTRUE},
 520       {kTRUE,kTRUE,kTRUE,kTRUE},
 521       {kTRUE,kTRUE,kTRUE,kTRUE},
 522       {kTRUE,kTRUE,kTRUE,kTRUE},
 523       {kTRUE,kTRUE,kTRUE,kTRUE},
 524       {kTRUE,kTRUE,kTRUE,kTRUE},
 525       {kTRUE,kTRUE,kTRUE,kTRUE},
 526       {kTRUE,kTRUE,kTRUE,kTRUE},
 527       {kTRUE,kTRUE,kTRUE,kTRUE},
 528       {kTRUE,kTRUE,kTRUE,kTRUE},
 529       {kFALSE,kFALSE,kFALSE,kFALSE},
 530       {kFALSE,kFALSE,kTRUE,kTRUE},
 531       {kTRUE,kTRUE,kTRUE,kTRUE},
 532       {kTRUE,kTRUE,kTRUE,kTRUE},
 533       {kTRUE,kTRUE,kFALSE,kFALSE},
 534       {kTRUE,kTRUE,kTRUE,kTRUE},
 535       {kTRUE,kTRUE,kTRUE,kTRUE},
 536       {kTRUE,kTRUE,kTRUE,kTRUE},
 537       {kFALSE,kFALSE,kFALSE,kFALSE},
 538       {kFALSE,kFALSE,kFALSE,kFALSE},
 539       {kTRUE,kTRUE,kTRUE,kTRUE},
 540       {kTRUE,kTRUE,kTRUE,kTRUE},
 541       {kFALSE,kFALSE,kFALSE,kFALSE},
 542       {kFALSE,kFALSE,kFALSE,kFALSE},
 543       {kTRUE,kTRUE,kTRUE,kTRUE},
 544       {kTRUE,kTRUE,kTRUE,kTRUE},
 545       {kFALSE,kFALSE,kFALSE,kFALSE},
 546       {kFALSE,kFALSE,kFALSE,kFALSE},
 547       {kFALSE,kFALSE,kFALSE,kFALSE},
 548       {kFALSE,kFALSE,kFALSE,kFALSE},
 549       {kTRUE,kTRUE,kTRUE,kTRUE},
 550       {kTRUE,kTRUE,kTRUE,kTRUE},
 551       {kTRUE,kTRUE,kFALSE,kFALSE},
 552       {kTRUE,kTRUE,kTRUE,kTRUE},
 553       {kFALSE,kFALSE,kFALSE,kFALSE},
 554       {kFALSE,kFALSE,kFALSE,kFALSE},
 555       {kFALSE,kFALSE,kFALSE,kFALSE},
 556       {kFALSE,kFALSE,kFALSE,kFALSE}
 557     };
 558     Double_t xyz1[3],xyz2[3];
 559     esdTrack.GetXYZAt(3.9,0.,xyz1);
 560     esdTrack.GetXYZAt(7.6,0.,xyz2);
 561     Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
 562     if(phi1<0) phi1+=2*TMath::Pi();
 563     Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
 564     Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
 565     if(phi2<0) phi2+=2*TMath::Pi();
 566     Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
 567     Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
 568     Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
 569     Bool_t lay1ok=kFALSE;
 570     if(mod1>=0 && mod1<4 && lad1<20){
 571       lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
 572     }
 573     Bool_t lay2ok=kFALSE;
 574     if(mod2>=0 && mod2<4 && lad2<40){
 575       lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
 576     }
 577     if(!lay1ok && !lay2ok) retval=kFALSE;
 578   }
 579
 580   return retval;
 581 }
 582 //---------------------------------------------------------------------------
 583 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
 584   // Set the pt bins
 585
 586   if(fPtBinLimits) {
 587     delete [] fPtBinLimits;
 588     fPtBinLimits = NULL;
 589     printf("Changing the pt bins\n");
 590   }
 591
 592   if(nPtBinLimits != fnPtBins+1){
 593     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
 594     SetNPtBins(nPtBinLimits-1);
 595   }
 596
 597   fnPtBinLimits = nPtBinLimits;
 598   SetGlobalIndex();
 599   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
 600   fPtBinLimits = new Float_t[fnPtBinLimits];
 601   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
 602
 603   return;
 604 }
 605 //---------------------------------------------------------------------------
 606 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
 607   // Set the variable names
 608
 609   if(fVarNames) {
 610     delete [] fVarNames;
 611     fVarNames = NULL;
 612     //printf("Changing the variable names\n");
 613   }
 614   if(nVars!=fnVars){
 615     printf("Wrong number of variables: it has to be %d\n",fnVars);
 616     return;
 617   }
 618   //fnVars=nVars;
 619   fVarNames = new TString[nVars];
 620   fIsUpperCut = new Bool_t[nVars];
 621   for(Int_t iv=0; iv<nVars; iv++) {
 622     fVarNames[iv] = varNames[iv];
 623     fIsUpperCut[iv] = isUpperCut[iv];
 624   }
 625
 626   return;
 627 }
 628 //---------------------------------------------------------------------------
 629 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
 630   // Set the variables to be used for cuts optimization
 631
 632   if(fVarsForOpt) {
 633     delete [] fVarsForOpt;
 634     fVarsForOpt = NULL;
 635     //printf("Changing the variables for cut optimization\n");
 636   }
 637
 638   if(nVars==0){//!=fnVars) {
 639     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
 640     return;
 641   }
 642
 643   fnVarsForOpt = 0;
 644   fVarsForOpt = new Bool_t[fnVars];
 645   for(Int_t iv=0; iv<fnVars; iv++) {
 646     fVarsForOpt[iv]=forOpt[iv];
 647     if(fVarsForOpt[iv]) fnVarsForOpt++;
 648   }
 649
 650   return;
 651 }
 652
 653 //---------------------------------------------------------------------------
 654 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
 655   //
 656   // set centrality estimator
 657   //
 658   fUseCentrality=flag;
 659   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
 660
 661   return;
 662 }
 663
 664
 665 //---------------------------------------------------------------------------
 666 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
 667   //
 668   // store the cuts
 669   //
 670   if(nVars!=fnVars) {
 671     printf("Wrong number of variables: it has to be %d\n",fnVars);
 672     AliFatal("exiting");
 673   }
 674   if(nPtBins!=fnPtBins) {
 675     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
 676     AliFatal("exiting");
 677   }
 678
 679   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 680
 681
 682   for(Int_t iv=0; iv<fnVars; iv++) {
 683
 684     for(Int_t ib=0; ib<fnPtBins; ib++) {
 685
 686       //check
 687       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
 688         cout<<"Overflow, exit..."<<endl;
 689         return;
 690       }
 691
 692       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
 693
 694     }
 695   }
 696   return;
 697 }
 698 //---------------------------------------------------------------------------
 699 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
 700   //
 701   // store the cuts
 702   //
 703   if(glIndex != fGlobalIndex){
 704     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
 705     AliFatal("exiting");
 706   }
 707   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 708
 709   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
 710     fCutsRD[iGl] = cutsRDGlob[iGl];
 711   }
 712   return;
 713 }
 714 //---------------------------------------------------------------------------
 715 void AliRDHFCuts::PrintAll() const {
 716   //
 717   // print all cuts values
 718   //
 719
 720   printf("Minimum vtx type %d\n",fMinVtxType);
 721   printf("Minimum vtx contr %d\n",fMinVtxContr);
 722   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
 723   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
 724   printf("Use PID %d\n",(Int_t)fUsePID);
 725   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
 726   printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
 727   printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
 728   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
 729   if(fOptPileup==1) printf(" -- Reject pileup event");
 730   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
 731   if(fUseCentrality>0) {
 732     TString estimator="";
 733     if(fUseCentrality==1) estimator = "V0";
 734     if(fUseCentrality==2) estimator = "Tracks";
 735     if(fUseCentrality==3) estimator = "Tracklets";
 736     if(fUseCentrality==4) estimator = "SPD clusters outer";
 737     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
 738   }
 739   if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
 740
 741   if(fVarNames){
 742     cout<<"Array of variables"<<endl;
 743     for(Int_t iv=0;iv<fnVars;iv++){
 744       cout<<fVarNames[iv]<<"\t";
 745     }
 746     cout<<endl;
 747   }
 748   if(fVarsForOpt){
 749     cout<<"Array of optimization"<<endl;
 750     for(Int_t iv=0;iv<fnVars;iv++){
 751       cout<<fVarsForOpt[iv]<<"\t";
 752     }
 753     cout<<endl;
 754   }
 755   if(fIsUpperCut){
 756     cout<<"Array of upper/lower cut"<<endl;
 757    for(Int_t iv=0;iv<fnVars;iv++){
 758      cout<<fIsUpperCut[iv]<<"\t";
 759    }
 760    cout<<endl;
 761   }
 762   if(fPtBinLimits){
 763     cout<<"Array of ptbin limits"<<endl;
 764     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
 765       cout<<fPtBinLimits[ib]<<"\t";
 766     }
 767     cout<<endl;
 768   }
 769   if(fCutsRD){
 770     cout<<"Matrix of cuts"<<endl;
 771    for(Int_t iv=0;iv<fnVars;iv++){
 772      for(Int_t ib=0;ib<fnPtBins;ib++){
 773        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
 774      }
 775      cout<<endl;
 776    }
 777    cout<<endl;
 778   }
 779   return;
 780 }
 781
 782 //--------------------------------------------------------------------------
 783 void AliRDHFCuts::PrintTrigger() const{
 784   // print the trigger selection
 785
 786   printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
 787
 788   cout<<" Trigger selection pattern: ";
 789   if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
 790   if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
 791   if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
 792   if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
 793   if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
 794   if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
 795   if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
 796   if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
 797   if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
 798   if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
 799   cout << endl<< endl;
 800
 801 }
 802
 803 //---------------------------------------------------------------------------
 804 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
 805   //
 806   // get the cuts
 807   //
 808
 809   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
 810
 811
 812   Int_t iv,ib;
 813   if(!cutsRD) {
 814     //cout<<"Initialization..."<<endl;
 815     cutsRD=new Float_t*[fnVars];
 816     for(iv=0; iv<fnVars; iv++) {
 817       cutsRD[iv] = new Float_t[fnPtBins];
 818     }
 819   }
 820
 821   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
 822     GetVarPtIndex(iGlobal,iv,ib);
 823     cutsRD[iv][ib] = fCutsRD[iGlobal];
 824   }
 825
 826   return;
 827 }
 828
 829 //---------------------------------------------------------------------------
 830 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
 831   //
 832   // give the global index from variable and pt bin
 833   //
 834   return iPtBin*fnVars+iVar;
 835 }
 836
 837 //---------------------------------------------------------------------------
 838 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
 839   //
 840   //give the index of the variable and of the pt bin from the global index
 841   //
 842   iPtBin=(Int_t)iGlob/fnVars;
 843   iVar=iGlob%fnVars;
 844
 845   return;
 846 }
 847
 848 //---------------------------------------------------------------------------
 849 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
 850   //
 851   //give the pt bin where the pt lies.
 852   //
 853   Int_t ptbin=-1;
 854   if(pt<fPtBinLimits[0])return ptbin;
 855   for (Int_t i=0;i<fnPtBins;i++){
 856     if(pt<fPtBinLimits[i+1]) {
 857       ptbin=i;
 858       break;
 859     }
 860   }
 861   return ptbin;
 862 }
 863 //-------------------------------------------------------------------
 864 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
 865   //
 866   // Give the value of cut set for the variable iVar and the pt bin iPtBin
 867   //
 868   if(!fCutsRD){
 869     cout<<"Cuts not iniziaisez yet"<<endl;
 870     return 0;
 871   }
 872   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
 873 }
 874 //-------------------------------------------------------------------
 875 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
 876   //
 877   // Get centrality percentile
 878   //
 879
 880   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 881   if(mcArray) {fUseAOD049=kFALSE;}
 882
 883   AliAODHeader *header=aodEvent->GetHeader();
 884   AliCentrality *centrality=header->GetCentralityP();
 885   Float_t cent=-999.;
 886   Bool_t isSelRun=kFALSE;
 887   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
 888   if(!centrality) return cent;
 889   else{
 890     if (estimator==kCentV0M){
 891       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
 892       if(cent<0){
 893         Int_t quality = centrality->GetQuality();
 894         if(quality<=1){
 895           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 896         }else{
 897           Int_t runnum=aodEvent->GetRunNumber();
 898           for(Int_t ir=0;ir<5;ir++){
 899             if(runnum==selRun[ir]){
 900               isSelRun=kTRUE;
 901               break;
 902             }
 903           }
 904           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 905         }
 906       }
 907
 908       //temporary fix for AOD049 outliers
 909       if(fUseAOD049&&cent>=0){
 910         Float_t v0=0;
 911         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
 912         v0+=aodV0->GetMTotV0A();
 913         v0+=aodV0->GetMTotV0C();
 914         if(cent==0&&v0<19500)return -1;//filtering issue
 915         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
 916         Float_t val= 1.30552 +  0.147931 * v0;
 917         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};
 918         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
 919       }
 920     }
 921     else {
 922       if (estimator==kCentTRK) {
 923         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
 924         if(cent<0){
 925           Int_t quality = centrality->GetQuality();
 926           if(quality<=1){
 927             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
 928           }else{
 929             Int_t runnum=aodEvent->GetRunNumber();
 930             for(Int_t ir=0;ir<5;ir++){
 931               if(runnum==selRun[ir]){
 932                 isSelRun=kTRUE;
 933                 break;
 934               }
 935             }
 936             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
 937           }
 938         }
 939       }
 940       else{
 941         if (estimator==kCentTKL){
 942           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
 943           if(cent<0){
 944             Int_t quality = centrality->GetQuality();
 945             if(quality<=1){
 946               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
 947             }else{
 948               Int_t runnum=aodEvent->GetRunNumber();
 949               for(Int_t ir=0;ir<5;ir++){
 950                 if(runnum==selRun[ir]){
 951                   isSelRun=kTRUE;
 952                   break;
 953             }
 954               }
 955               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
 956             }
 957           }
 958         }
 959         else{
 960           if (estimator==kCentCL1){
 961             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
 962             if(cent<0){
 963               Int_t quality = centrality->GetQuality();
 964               if(quality<=1){
 965                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
 966               }else{
 967                 Int_t runnum=aodEvent->GetRunNumber();
 968                 for(Int_t ir=0;ir<5;ir++){
 969                   if(runnum==selRun[ir]){
 970                     isSelRun=kTRUE;
 971                     break;
 972                   }
 973                 }
 974                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
 975               }
 976             }
 977           }
 978           else {
 979             AliWarning("Centrality estimator not valid");
 980
 981           }
 982         }
 983       }
 984     }
 985   }
 986   return cent;
 987 }
 988 //-------------------------------------------------------------------
 989 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
 990   //
 991   // Compare two cuts objects
 992   //
 993
 994   Bool_t areEqual=kTRUE;
 995
 996   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
 997
 998   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
 999
1000   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1001
1002   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1003
1004   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1005
1006   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1007   if(fTrackCuts){
1008     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1009
1010     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1011
1012     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1013
1014     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;}
1015   }
1016
1017   if(fCutsRD) {
1018    for(Int_t iv=0;iv<fnVars;iv++) {
1019      for(Int_t ib=0;ib<fnPtBins;ib++) {
1020        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1021          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1022          areEqual=kFALSE;
1023        }
1024      }
1025    }
1026   }
1027
1028   return areEqual;
1029 }
1030 //---------------------------------------------------------------------------
1031 void AliRDHFCuts::MakeTable() const {
1032   //
1033   // print cuts values in table format
1034   //
1035
1036         TString ptString = "pT range";
1037         if(fVarNames && fPtBinLimits && fCutsRD){
1038                 TString firstLine(Form("*       %-15s",ptString.Data()));
1039                 for (Int_t ivar=0; ivar<fnVars; ivar++){
1040                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
1041                         if (ivar == fnVars){
1042                                 firstLine+="*\n";
1043                         }
1044                 }
1045                 Printf("%s",firstLine.Data());
1046
1047                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1048                         TString line;
1049                         if (ipt==fnPtBins-1){
1050                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
1051                         }
1052                         else{
1053                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1054                         }
1055                         for (Int_t ivar=0; ivar<fnVars; ivar++){
1056                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1057                         }
1058                         Printf("%s",line.Data());
1059                 }
1060
1061         }
1062
1063
1064   return;
1065 }
1066 //--------------------------------------------------------------------------
1067 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1068                                         AliAODEvent *aod) const
1069 {
1070   //
1071   // Recalculate primary vertex without daughters
1072   //
1073
1074   if(!aod) {
1075     AliError("Can not remove daughters from vertex without AOD event");
1076     return 0;
1077   }
1078
1079   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1080   if(!recvtx){
1081     AliDebug(2,"Removal of daughter tracks failed");
1082     return kFALSE;
1083   }
1084
1085
1086   //set recalculed primary vertex
1087   d->SetOwnPrimaryVtx(recvtx);
1088   delete recvtx;
1089
1090   return kTRUE;
1091 }
1092 //--------------------------------------------------------------------------
1093 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1094 {
1095   //
1096   // Recalculate primary vertex without daughters
1097   //
1098
1099   if(!aod) {
1100     AliError("Can not get MC vertex without AOD event");
1101     return kFALSE;
1102   }
1103
1104   // load MC header
1105   AliAODMCHeader *mcHeader =
1106     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1107   if(!mcHeader) {
1108     AliError("Can not get MC vertex without AODMCHeader event");
1109     return kFALSE;
1110   }
1111   Double_t pos[3];
1112   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1113   mcHeader->GetVertex(pos);
1114   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1115
1116   if(!recvtx){
1117     AliDebug(2,"Removal of daughter tracks failed");
1118     return kFALSE;
1119   }
1120
1121   //set recalculed primary vertex
1122   d->SetOwnPrimaryVtx(recvtx);
1123
1124   d->RecalculateImpPars(recvtx,aod);
1125
1126   delete recvtx;
1127
1128   return kTRUE;
1129 }
1130 //--------------------------------------------------------------------------
1131 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1132                                      AliAODEvent *aod,
1133                                      AliAODVertex *origownvtx) const
1134 {
1135   //
1136   // Clean-up own primary vertex if needed
1137   //
1138
1139   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1140     d->UnsetOwnPrimaryVtx();
1141     if(origownvtx) {
1142       d->SetOwnPrimaryVtx(origownvtx);
1143       delete origownvtx; origownvtx=NULL;
1144     }
1145     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1146   } else {
1147     if(origownvtx) {
1148       delete origownvtx; origownvtx=NULL;
1149     }
1150   }
1151   return;
1152 }
1153 //--------------------------------------------------------------------------
1154 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1155 {
1156   //
1157   // Checks if this candidate is matched to MC signal
1158   //
1159
1160   if(!aod) return kFALSE;
1161
1162   // get the MC array
1163   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1164
1165   if(!mcArray) return kFALSE;
1166
1167   // try to match
1168   Int_t label = d->MatchToMC(pdg,mcArray);
1169
1170   if(label>=0) {
1171     //printf("MATCH!\n");
1172     return kTRUE;
1173   }
1174
1175   return kFALSE;
1176 }
1177
1178
1179 //--------------------------------------------------------------------------
1180 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1181   // recompute event primary vertex from AOD tracks
1182
1183    AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1184    vertexer->SetITSMode();
1185    vertexer->SetMinClusters(3);
1186
1187    AliAODVertex* pvtx=event->GetPrimaryVertex();
1188    if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1189      Float_t diamondcovxy[3];
1190      event->GetDiamondCovXY(diamondcovxy);
1191      Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1192      Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1193      AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1194      vertexer->SetVtxStart(diamond);
1195      delete diamond; diamond=NULL;
1196    }
1197
1198    AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1199    if(!vertexESD) return kFALSE;
1200    if(vertexESD->GetNContributors()<=0) {
1201      //AliDebug(2,"vertexing failed");
1202      delete vertexESD; vertexESD=NULL;
1203      return kFALSE;
1204    }
1205    delete vertexer; vertexer=NULL;
1206
1207    // convert to AliAODVertex
1208    Double_t pos[3],cov[6],chi2perNDF;
1209    vertexESD->GetXYZ(pos); // position
1210    vertexESD->GetCovMatrix(cov); //covariance matrix
1211    chi2perNDF = vertexESD->GetChi2toNDF();
1212    delete vertexESD; vertexESD=NULL;
1213
1214    pvtx->SetPosition(pos[0],pos[1],pos[2]);
1215    pvtx->SetChi2perNDF(chi2perNDF);
1216    pvtx->SetCovMatrix(cov);
1217
1218    return kTRUE;
1219 }