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