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