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