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