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