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