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