PWG3/vertexingHF/AliRDHFCuts.cxx

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