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