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