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