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