PWG3/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 "AliRDHFCuts.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliInputEventHandler.h"
  44 #include "AliPIDResponse.h"
  45
  46 ClassImp(AliRDHFCuts)
  47
  48 //--------------------------------------------------------------------------
  49 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
  50 AliAnalysisCuts(name,title),
  51 fMinVtxType(3),
  52 fMinVtxContr(1),
  53 fMaxVtxRedChi2(1e6),
  54 fMaxVtxZ(10.),
  55 fMinSPDMultiplicity(0),
  56 fTriggerMask(0),
  57 fTriggerClass("CINT1"),
  58 fTrackCuts(0),
  59 fnPtBins(1),
  60 fnPtBinLimits(1),
  61 fPtBinLimits(0),
  62 fnVars(1),
  63 fVarNames(0),
  64 fnVarsForOpt(0),
  65 fVarsForOpt(0),
  66 fGlobalIndex(1),
  67 fCutsRD(0),
  68 fIsUpperCut(0),
  69 fUsePID(kFALSE),
  70 fUseAOD049(kFALSE),
  71 fPidHF(0),
  72 fWhyRejection(0),
  73 fRemoveDaughtersFromPrimary(kFALSE),
  74 fUseMCVertex(kFALSE),
  75 fOptPileup(0),
  76 fMinContrPileup(3),
  77 fMinDzPileup(0.6),
  78 fUseCentrality(0),
  79 fMinCentrality(0.),
  80 fMaxCentrality(100.),
  81 fFixRefs(kFALSE),
  82 fIsSelectedCuts(0),
  83 fIsSelectedPID(0),
  84 fMinPtCand(-1.),
  85 fMaxPtCand(100000.)
  86 {
  87   //
  88   // Default Constructor
  89   //
  90 }
  91 //--------------------------------------------------------------------------
  92 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
  93   AliAnalysisCuts(source),
  94   fMinVtxType(source.fMinVtxType),
  95   fMinVtxContr(source.fMinVtxContr),
  96   fMaxVtxRedChi2(source.fMaxVtxRedChi2),
  97   fMaxVtxZ(source.fMaxVtxZ),
  98   fMinSPDMultiplicity(source.fMinSPDMultiplicity),
  99   fTriggerMask(source.fTriggerMask),
 100   fTriggerClass(source.fTriggerClass),
 101   fTrackCuts(0),
 102   fnPtBins(source.fnPtBins),
 103   fnPtBinLimits(source.fnPtBinLimits),
 104   fPtBinLimits(0),
 105   fnVars(source.fnVars),
 106   fVarNames(0),
 107   fnVarsForOpt(source.fnVarsForOpt),
 108   fVarsForOpt(0),
 109   fGlobalIndex(source.fGlobalIndex),
 110   fCutsRD(0),
 111   fIsUpperCut(0),
 112   fUsePID(source.fUsePID),
 113   fUseAOD049(source.fUseAOD049),
 114   fPidHF(0),
 115   fWhyRejection(source.fWhyRejection),
 116   fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
 117   fUseMCVertex(source.fUseMCVertex),
 118   fOptPileup(source.fOptPileup),
 119   fMinContrPileup(source.fMinContrPileup),
 120   fMinDzPileup(source.fMinDzPileup),
 121   fUseCentrality(source.fUseCentrality),
 122   fMinCentrality(source.fMinCentrality),
 123   fMaxCentrality(source.fMaxCentrality),
 124   fFixRefs(source.fFixRefs),
 125   fIsSelectedCuts(source.fIsSelectedCuts),
 126   fIsSelectedPID(source.fIsSelectedPID),
 127   fMinPtCand(source.fMinPtCand),
 128   fMaxPtCand(source.fMaxPtCand)
 129 {
 130   //
 131   // Copy constructor
 132   //
 133   cout<<"Copy constructor"<<endl;
 134   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
 135   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 136   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 137   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 138   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 139   if(source.fPidHF) SetPidHF(source.fPidHF);
 140   PrintAll();
 141
 142 }
 143 //--------------------------------------------------------------------------
 144 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
 145 {
 146   //
 147   // assignment operator
 148   //
 149   if(&source == this) return *this;
 150
 151   AliAnalysisCuts::operator=(source);
 152
 153   fMinVtxType=source.fMinVtxType;
 154   fMinVtxContr=source.fMinVtxContr;
 155   fMaxVtxRedChi2=source.fMaxVtxRedChi2;
 156   fMaxVtxZ=source.fMaxVtxZ;
 157   fMinSPDMultiplicity=source.fMinSPDMultiplicity;
 158   fTriggerMask=source.fTriggerMask;
 159   fTriggerClass=source.fTriggerClass;
 160   fnPtBins=source.fnPtBins;
 161   fnVars=source.fnVars;
 162   fGlobalIndex=source.fGlobalIndex;
 163   fnVarsForOpt=source.fnVarsForOpt;
 164   fUsePID=source.fUsePID;
 165   fUseAOD049=source.fUseAOD049;
 166   SetPidHF(source.GetPidHF());
 167   fWhyRejection=source.fWhyRejection;
 168   fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
 169   fUseMCVertex=source.fUseMCVertex;
 170   fOptPileup=source.fOptPileup;
 171   fMinContrPileup=source.fMinContrPileup;
 172   fMinDzPileup=source.fMinDzPileup;
 173   fUseCentrality=source.fUseCentrality;
 174   fMinCentrality=source.fMinCentrality;
 175   fMaxCentrality=source.fMaxCentrality;
 176   fFixRefs=source.fFixRefs;
 177   fIsSelectedCuts=source.fIsSelectedCuts;
 178   fIsSelectedPID=source.fIsSelectedPID;
 179   fMinPtCand=source.fMinPtCand;
 180   fMaxPtCand=source.fMaxPtCand;
 181
 182   if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
 183   if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
 184   if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
 185   if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
 186   if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
 187   PrintAll();
 188
 189   return *this;
 190 }
 191 //--------------------------------------------------------------------------
 192 AliRDHFCuts::~AliRDHFCuts() {
 193   //
 194   // Default Destructor
 195   //
 196   if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
 197   if(fVarNames) {delete [] fVarNames; fVarNames=0;}
 198   if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
 199   if(fCutsRD) {
 200     delete [] fCutsRD;
 201     fCutsRD=0;
 202   }
 203   if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
 204   if(fPidHF){
 205     delete fPidHF;
 206     fPidHF=0;
 207   }
 208 }
 209 //---------------------------------------------------------------------------
 210 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
 211   //
 212   // Centrality selection
 213   //
 214
 215   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
 216     AliWarning("Centrality estimator not valid");
 217     return 3;
 218   }else{
 219     Float_t centvalue=GetCentrality((AliAODEvent*)event);
 220     if (centvalue<-998.){//-999 if no centralityP
 221       return 0;
 222     }else{
 223       if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
 224         return 2;
 225       }
 226     }
 227   }
 228   return 0;
 229 }
 230 //---------------------------------------------------------------------------
 231 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
 232   //
 233   // Event selection
 234   //
 235   //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
 236
 237   fWhyRejection=0;
 238
 239   // check if it's MC
 240   Bool_t isMC=kFALSE;
 241   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 242   if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
 243
 244   // settings for the TPC dE/dx BB parameterization
 245   if(fPidHF && fPidHF->GetOldPid()) {
 246     // pp, from LHC10d onwards
 247     if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
 248        event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
 249     // PbPb LHC10h
 250     if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
 251     // MC
 252     if(isMC) fPidHF->SetMC(kTRUE);
 253     if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
 254       fPidHF->SetMClowenpp2011(kTRUE);
 255   }
 256   else if(fPidHF && !fPidHF->GetOldPid()) {
 257    if(fPidHF->GetPidResponse()==0x0){
 258     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 259     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 260     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 261     fPidHF->SetPidResponse(pidResp);
 262    }
 263   }
 264
 265
 266   // trigger class
 267   TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
 268   // don't do for MC and for PbPb 2010 data
 269   if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
 270     if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {
 271       fWhyRejection=5;
 272       return kFALSE;
 273     }
 274   }
 275
 276   // TEMPORARY FIX FOR REFERENCES
 277   // Fix references to daughter tracks
 278   //  if(fFixRefs) {
 279   //    AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
 280   //    fixer->FixReferences((AliAODEvent*)event);
 281   //    delete fixer;
 282   //  }
 283   //
 284
 285
 286
 287   // vertex requirements
 288
 289   const AliVVertex *vertex = event->GetPrimaryVertex();
 290
 291   if(!vertex) return kFALSE;
 292
 293   TString title=vertex->GetTitle();
 294   if(title.Contains("Z") && fMinVtxType>1) return kFALSE;
 295   if(title.Contains("3D") && fMinVtxType>2) return kFALSE;
 296
 297   if(vertex->GetNContributors()<fMinVtxContr) return kFALSE;
 298
 299   if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
 300     fWhyRejection=6;
 301     return kFALSE;
 302   }
 303
 304
 305   // pile-up rejection
 306   if(fOptPileup==kRejectPileupEvent){
 307     Int_t cutc=(Int_t)fMinContrPileup;
 308     Double_t cutz=(Double_t)fMinDzPileup;
 309     if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
 310       fWhyRejection=1;
 311       return kFALSE;
 312     }
 313   }
 314
 315   // centrality selection
 316   if (fUseCentrality!=kCentOff) {
 317     Int_t rejection=IsEventSelectedInCentrality(event);
 318     if(rejection>1){
 319       fWhyRejection=rejection;
 320       return kFALSE;
 321     }
 322   }
 323
 324   return kTRUE;
 325 }
 326 //---------------------------------------------------------------------------
 327 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
 328   //
 329   // Daughter tracks selection
 330   //
 331   if(!fTrackCuts) return kTRUE;
 332
 333   Int_t ndaughters = d->GetNDaughters();
 334   AliAODVertex *vAOD = d->GetPrimaryVtx();
 335   Double_t pos[3],cov[6];
 336   vAOD->GetXYZ(pos);
 337   vAOD->GetCovarianceMatrix(cov);
 338   const AliESDVertex vESD(pos,cov,100.,100);
 339
 340   Bool_t retval=kTRUE;
 341
 342   for(Int_t idg=0; idg<ndaughters; idg++) {
 343     AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
 344     if(!dgTrack) {retval = kFALSE; continue;}
 345     //printf("charge %d\n",dgTrack->Charge());
 346     if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
 347
 348     if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
 349   }
 350
 351   return retval;
 352 }
 353 //---------------------------------------------------------------------------
 354 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
 355   //
 356   // Convert to ESDtrack, relate to vertex and check cuts
 357   //
 358   if(!cuts) return kTRUE;
 359
 360   Bool_t retval=kTRUE;
 361
 362   // convert to ESD track here
 363   AliESDtrack esdTrack(track);
 364   // set the TPC cluster info
 365   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
 366   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
 367   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
 368   // needed to calculate the impact parameters
 369   esdTrack.RelateToVertex(primary,0.,3.);
 370
 371   if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
 372
 373   if(fOptPileup==kRejectTracksFromPileupVertex){
 374     // to be implemented
 375     // we need either to have here the AOD Event,
 376     // or to have the pileup vertex object
 377   }
 378   return retval;
 379 }
 380 //---------------------------------------------------------------------------
 381 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
 382   // Set the pt bins
 383
 384   if(fPtBinLimits) {
 385     delete [] fPtBinLimits;
 386     fPtBinLimits = NULL;
 387     printf("Changing the pt bins\n");
 388   }
 389
 390   if(nPtBinLimits != fnPtBins+1){
 391     cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
 392     SetNPtBins(nPtBinLimits-1);
 393   }
 394
 395   fnPtBinLimits = nPtBinLimits;
 396   SetGlobalIndex();
 397   //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
 398   fPtBinLimits = new Float_t[fnPtBinLimits];
 399   for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
 400
 401   return;
 402 }
 403 //---------------------------------------------------------------------------
 404 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
 405   // Set the variable names
 406
 407   if(fVarNames) {
 408     delete [] fVarNames;
 409     fVarNames = NULL;
 410     //printf("Changing the variable names\n");
 411   }
 412   if(nVars!=fnVars){
 413     printf("Wrong number of variables: it has to be %d\n",fnVars);
 414     return;
 415   }
 416   //fnVars=nVars;
 417   fVarNames = new TString[nVars];
 418   fIsUpperCut = new Bool_t[nVars];
 419   for(Int_t iv=0; iv<nVars; iv++) {
 420     fVarNames[iv] = varNames[iv];
 421     fIsUpperCut[iv] = isUpperCut[iv];
 422   }
 423
 424   return;
 425 }
 426 //---------------------------------------------------------------------------
 427 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
 428   // Set the variables to be used for cuts optimization
 429
 430   if(fVarsForOpt) {
 431     delete [] fVarsForOpt;
 432     fVarsForOpt = NULL;
 433     //printf("Changing the variables for cut optimization\n");
 434   }
 435
 436   if(nVars==0){//!=fnVars) {
 437     printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
 438     return;
 439   }
 440
 441   fnVarsForOpt = 0;
 442   fVarsForOpt = new Bool_t[fnVars];
 443   for(Int_t iv=0; iv<fnVars; iv++) {
 444     fVarsForOpt[iv]=forOpt[iv];
 445     if(fVarsForOpt[iv]) fnVarsForOpt++;
 446   }
 447
 448   return;
 449 }
 450
 451 //---------------------------------------------------------------------------
 452 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
 453   //
 454   // set centrality estimator
 455   //
 456   fUseCentrality=flag;
 457   if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
 458
 459   return;
 460 }
 461
 462
 463 //---------------------------------------------------------------------------
 464 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
 465   //
 466   // store the cuts
 467   //
 468   if(nVars!=fnVars) {
 469     printf("Wrong number of variables: it has to be %d\n",fnVars);
 470     return;
 471   }
 472   if(nPtBins!=fnPtBins) {
 473     printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
 474     return;
 475   }
 476
 477   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 478
 479
 480   for(Int_t iv=0; iv<fnVars; iv++) {
 481
 482     for(Int_t ib=0; ib<fnPtBins; ib++) {
 483
 484       //check
 485       if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
 486         cout<<"Overflow, exit..."<<endl;
 487         return;
 488       }
 489
 490       fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
 491
 492     }
 493   }
 494   return;
 495 }
 496 //---------------------------------------------------------------------------
 497 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
 498   //
 499   // store the cuts
 500   //
 501   if(glIndex != fGlobalIndex){
 502     cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
 503     return;
 504   }
 505   if(!fCutsRD)  fCutsRD = new Float_t[fGlobalIndex];
 506
 507   for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
 508     fCutsRD[iGl] = cutsRDGlob[iGl];
 509   }
 510   return;
 511 }
 512 //---------------------------------------------------------------------------
 513 void AliRDHFCuts::PrintAll() const {
 514   //
 515   // print all cuts values
 516   //
 517
 518   printf("Minimum vtx type %d\n",fMinVtxType);
 519   printf("Minimum vtx contr %d\n",fMinVtxContr);
 520   printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
 521   printf("Min SPD mult %d\n",fMinSPDMultiplicity);
 522   printf("Use PID %d\n",(Int_t)fUsePID);
 523   printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
 524   printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
 525   if(fOptPileup==1) printf(" -- Reject pileup event");
 526   if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
 527   if(fUseCentrality>0) {
 528     TString estimator="";
 529     if(fUseCentrality==1) estimator = "V0";
 530     if(fUseCentrality==2) estimator = "Tracks";
 531     if(fUseCentrality==3) estimator = "Tracklets";
 532     if(fUseCentrality==4) estimator = "SPD clusters outer";
 533     printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
 534   }
 535
 536   if(fVarNames){
 537     cout<<"Array of variables"<<endl;
 538     for(Int_t iv=0;iv<fnVars;iv++){
 539       cout<<fVarNames[iv]<<"\t";
 540     }
 541     cout<<endl;
 542   }
 543   if(fVarsForOpt){
 544     cout<<"Array of optimization"<<endl;
 545     for(Int_t iv=0;iv<fnVars;iv++){
 546       cout<<fVarsForOpt[iv]<<"\t";
 547     }
 548     cout<<endl;
 549   }
 550   if(fIsUpperCut){
 551     cout<<"Array of upper/lower cut"<<endl;
 552    for(Int_t iv=0;iv<fnVars;iv++){
 553      cout<<fIsUpperCut[iv]<<"\t";
 554    }
 555    cout<<endl;
 556   }
 557   if(fPtBinLimits){
 558     cout<<"Array of ptbin limits"<<endl;
 559     for(Int_t ib=0;ib<fnPtBinLimits;ib++){
 560       cout<<fPtBinLimits[ib]<<"\t";
 561     }
 562     cout<<endl;
 563   }
 564   if(fCutsRD){
 565     cout<<"Matrix of cuts"<<endl;
 566    for(Int_t iv=0;iv<fnVars;iv++){
 567      for(Int_t ib=0;ib<fnPtBins;ib++){
 568        cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
 569      }
 570      cout<<endl;
 571    }
 572    cout<<endl;
 573   }
 574   return;
 575 }
 576 //---------------------------------------------------------------------------
 577 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
 578   //
 579   // get the cuts
 580   //
 581
 582   //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
 583
 584
 585   Int_t iv,ib;
 586   if(!cutsRD) {
 587     //cout<<"Initialization..."<<endl;
 588     cutsRD=new Float_t*[fnVars];
 589     for(iv=0; iv<fnVars; iv++) {
 590       cutsRD[iv] = new Float_t[fnPtBins];
 591     }
 592   }
 593
 594   for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
 595     GetVarPtIndex(iGlobal,iv,ib);
 596     cutsRD[iv][ib] = fCutsRD[iGlobal];
 597   }
 598
 599   return;
 600 }
 601
 602 //---------------------------------------------------------------------------
 603 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
 604   //
 605   // give the global index from variable and pt bin
 606   //
 607   return iPtBin*fnVars+iVar;
 608 }
 609
 610 //---------------------------------------------------------------------------
 611 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
 612   //
 613   //give the index of the variable and of the pt bin from the global index
 614   //
 615   iPtBin=(Int_t)iGlob/fnVars;
 616   iVar=iGlob%fnVars;
 617
 618   return;
 619 }
 620
 621 //---------------------------------------------------------------------------
 622 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
 623   //
 624   //give the pt bin where the pt lies.
 625   //
 626   Int_t ptbin=-1;
 627   if(pt<fPtBinLimits[0])return ptbin;
 628   for (Int_t i=0;i<fnPtBins;i++){
 629     if(pt<fPtBinLimits[i+1]) {
 630       ptbin=i;
 631       break;
 632     }
 633   }
 634   return ptbin;
 635 }
 636 //-------------------------------------------------------------------
 637 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
 638   //
 639   // Give the value of cut set for the variable iVar and the pt bin iPtBin
 640   //
 641   if(!fCutsRD){
 642     cout<<"Cuts not iniziaisez yet"<<endl;
 643     return 0;
 644   }
 645   return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
 646 }
 647 //-------------------------------------------------------------------
 648 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
 649   //
 650   // Get centrality percentile
 651   //
 652
 653   TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 654   if(mcArray) {fUseAOD049=kFALSE;}
 655
 656   AliAODHeader *header=aodEvent->GetHeader();
 657   AliCentrality *centrality=header->GetCentralityP();
 658   Float_t cent=-999.;
 659   Bool_t isSelRun=kFALSE;
 660   Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
 661   if(!centrality) return cent;
 662   else{
 663     if (estimator==kCentV0M){
 664       cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
 665       if(cent<0){
 666         Int_t quality = centrality->GetQuality();
 667         if(quality<=1){
 668           cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 669         }else{
 670           Int_t runnum=aodEvent->GetRunNumber();
 671           for(Int_t ir=0;ir<5;ir++){
 672             if(runnum==selRun[ir]){
 673               isSelRun=kTRUE;
 674               break;
 675             }
 676           }
 677           if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 678         }
 679       }
 680
 681       //temporary fix for AOD049 outliers
 682       if(fUseAOD049&&cent>=0){
 683         Float_t v0=0;
 684         AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
 685         v0+=aodV0->GetMTotV0A();
 686         v0+=aodV0->GetMTotV0C();
 687         if(cent==0&&v0<19500)return -1;//filtering issue
 688         Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
 689         Float_t val= 1.30552 +  0.147931 * v0;
 690         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};
 691         if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
 692       }
 693     }
 694     else {
 695       if (estimator==kCentTRK) {
 696         cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
 697         if(cent<0){
 698           Int_t quality = centrality->GetQuality();
 699           if(quality<=1){
 700             cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
 701           }else{
 702             Int_t runnum=aodEvent->GetRunNumber();
 703             for(Int_t ir=0;ir<5;ir++){
 704               if(runnum==selRun[ir]){
 705                 isSelRun=kTRUE;
 706                 break;
 707               }
 708             }
 709             if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
 710           }
 711         }
 712       }
 713       else{
 714         if (estimator==kCentTKL){
 715           cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
 716           if(cent<0){
 717             Int_t quality = centrality->GetQuality();
 718             if(quality<=1){
 719               cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
 720             }else{
 721               Int_t runnum=aodEvent->GetRunNumber();
 722               for(Int_t ir=0;ir<5;ir++){
 723                 if(runnum==selRun[ir]){
 724                   isSelRun=kTRUE;
 725                   break;
 726             }
 727               }
 728               if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
 729             }
 730           }
 731         }
 732         else{
 733           if (estimator==kCentCL1){
 734             cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
 735             if(cent<0){
 736               Int_t quality = centrality->GetQuality();
 737               if(quality<=1){
 738                 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
 739               }else{
 740                 Int_t runnum=aodEvent->GetRunNumber();
 741                 for(Int_t ir=0;ir<5;ir++){
 742                   if(runnum==selRun[ir]){
 743                     isSelRun=kTRUE;
 744                     break;
 745                   }
 746                 }
 747                 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
 748               }
 749             }
 750           }
 751           else {
 752             AliWarning("Centrality estimator not valid");
 753
 754           }
 755         }
 756       }
 757     }
 758   }
 759   return cent;
 760 }
 761 //-------------------------------------------------------------------
 762 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
 763   //
 764   // Compare two cuts objects
 765   //
 766
 767   Bool_t areEqual=kTRUE;
 768
 769   if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d  %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
 770
 771   if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d  %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
 772
 773   if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) {   printf("Max vtx red chi2 %f  %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
 774
 775   if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) {  printf("Min SPD mult %d\n  %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
 776
 777   if(fUsePID!=obj->fUsePID) { printf("Use PID %d  %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
 778
 779   if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d  %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
 780   if(fTrackCuts){
 781     if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d  %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
 782
 783     if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d  %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
 784
 785     if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f  %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
 786
 787     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;}
 788   }
 789
 790   if(fCutsRD) {
 791    for(Int_t iv=0;iv<fnVars;iv++) {
 792      for(Int_t ib=0;ib<fnPtBins;ib++) {
 793        if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
 794          cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"   "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
 795          areEqual=kFALSE;
 796        }
 797      }
 798    }
 799   }
 800
 801   return areEqual;
 802 }
 803 //---------------------------------------------------------------------------
 804 void AliRDHFCuts::MakeTable() const {
 805   //
 806   // print cuts values in table format
 807   //
 808
 809         TString ptString = "pT range";
 810         if(fVarNames && fPtBinLimits && fCutsRD){
 811                 TString firstLine(Form("*       %-15s",ptString.Data()));
 812                 for (Int_t ivar=0; ivar<fnVars; ivar++){
 813                         firstLine+=Form("*    %-15s  ",fVarNames[ivar].Data());
 814                         if (ivar == fnVars){
 815                                 firstLine+="*\n";
 816                         }
 817                 }
 818                 Printf("%s",firstLine.Data());
 819
 820                 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
 821                         TString line;
 822                         if (ipt==fnPtBins-1){
 823                                 line=Form("*  %5.1f < pt < inf    ",fPtBinLimits[ipt]);
 824                         }
 825                         else{
 826                                 line=Form("*  %5.1f < pt < %4.1f   ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
 827                         }
 828                         for (Int_t ivar=0; ivar<fnVars; ivar++){
 829                                 line+=Form("*     %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
 830                         }
 831                         Printf("%s",line.Data());
 832                 }
 833
 834         }
 835
 836
 837   return;
 838 }
 839 //--------------------------------------------------------------------------
 840 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
 841                                         AliAODEvent *aod) const
 842 {
 843   //
 844   // Recalculate primary vertex without daughters
 845   //
 846
 847   if(!aod) {
 848     AliError("Can not remove daughters from vertex without AOD event");
 849     return 0;
 850   }
 851
 852   AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
 853   if(!recvtx){
 854     AliDebug(2,"Removal of daughter tracks failed");
 855     return kFALSE;
 856   }
 857
 858
 859   //set recalculed primary vertex
 860   d->SetOwnPrimaryVtx(recvtx);
 861   delete recvtx;
 862
 863   return kTRUE;
 864 }
 865 //--------------------------------------------------------------------------
 866 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
 867 {
 868   //
 869   // Recalculate primary vertex without daughters
 870   //
 871
 872   if(!aod) {
 873     AliError("Can not get MC vertex without AOD event");
 874     return kFALSE;
 875   }
 876
 877   // load MC header
 878   AliAODMCHeader *mcHeader =
 879     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 880   if(!mcHeader) {
 881     AliError("Can not get MC vertex without AODMCHeader event");
 882     return kFALSE;
 883   }
 884   Double_t pos[3];
 885   Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
 886   mcHeader->GetVertex(pos);
 887   AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
 888
 889   if(!recvtx){
 890     AliDebug(2,"Removal of daughter tracks failed");
 891     return kFALSE;
 892   }
 893
 894   //set recalculed primary vertex
 895   d->SetOwnPrimaryVtx(recvtx);
 896
 897   d->RecalculateImpPars(recvtx,aod);
 898
 899   delete recvtx;
 900
 901   return kTRUE;
 902 }
 903 //--------------------------------------------------------------------------
 904 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
 905                                      AliAODEvent *aod,
 906                                      AliAODVertex *origownvtx) const
 907 {
 908   //
 909   // Clean-up own primary vertex if needed
 910   //
 911
 912   if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
 913     d->UnsetOwnPrimaryVtx();
 914     if(origownvtx) {
 915       d->SetOwnPrimaryVtx(origownvtx);
 916       delete origownvtx; origownvtx=NULL;
 917     }
 918     d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
 919   } else {
 920     if(origownvtx) {
 921       delete origownvtx; origownvtx=NULL;
 922     }
 923   }
 924   return;
 925 }
 926
 927
 928