AliPhysics  ff07904 (ff07904)
AliRDHFCuts.cxx
Go to the documentation of this file.
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 
19 //
20 // Base class for cuts on AOD reconstructed heavy-flavour decay
21 //
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
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 "AliAODHandler.h"
45 #include "AliInputEventHandler.h"
46 #include "AliPIDResponse.h"
47 #include "AliAnalysisUtils.h"
48 #include "AliMultSelection.h"
49 #include "AliAODVZERO.h"
50 #include "TRandom.h"
51 #include <TF1.h>
52 #include <TFile.h>
53 #include <TKey.h>
54 
55 using std::cout;
56 using std::endl;
57 
59 ClassImp(AliRDHFCuts);
61 
62 
63 //--------------------------------------------------------------------------
65 AliAnalysisCuts(name,title),
66 fMinVtxType(3),
67 fMinVtxContr(1),
68 fMaxVtxRedChi2(1e6),
69 fMaxVtxZ(10.),
70 fMinSPDMultiplicity(0),
71 fMaxVtxChi2PileupMV(5.),
72 fMinWDzPileupMV(15.),
73 fRejectPlpFromDiffBCMV(kFALSE),
74 fTriggerMask(AliVEvent::kAnyINT),
75 fUseOnlyOneTrigger(kFALSE),
76 fTrackCuts(0),
77 fnPtBins(1),
78 fnPtBinLimits(1),
79 fPtBinLimits(0),
80 fnVars(1),
81 fVarNames(0),
82 fnVarsForOpt(0),
83 fVarsForOpt(0),
84 fGlobalIndex(1),
85 fCutsRD(0),
86 fIsUpperCut(0),
87 fUsePID(kFALSE),
88 fUseAOD049(kFALSE),
89 fPidHF(0),
90 fWhyRejection(0),
91 fEvRejectionBits(0),
92 fRemoveDaughtersFromPrimary(kFALSE),
93 fUseMCVertex(kFALSE),
94 fUsePhysicsSelection(kTRUE),
95 fOptPileup(0),
96 fMinContrPileup(3),
97 fMinDzPileup(0.6),
98 fUseMultDepPileupCut(kFALSE),
99 fUseCentrality(0),
100 fMinCentrality(0.),
101 fMaxCentrality(100.),
102 fMultSelectionObjectName("MultSelection"),
103 fFixRefs(kFALSE),
104 fIsSelectedCuts(0),
105 fIsSelectedPID(0),
106 fMinPtCand(-1.),
107 fMaxPtCand(100000.),
108 fMaxRapidityCand(-999.),
109 fKeepSignalMC(kFALSE),
110 fIsCandTrackSPDFirst(kFALSE),
111 fMaxPtCandTrackSPDFirst(0.),
112 fApplySPDDeadPbPb2011(kFALSE),
113 fApplySPDMisalignedPP2012(kFALSE),
114 fMaxDiffTRKV0Centr(-1.),
115 fRemoveTrackletOutliers(kFALSE),
116 fCutOnzVertexSPD(3),
117 fKinkReject(kFALSE),
118 fUseTrackSelectionWithFilterBits(kTRUE),
119 fUseCentrFlatteningInMC(kFALSE),
120 fHistCentrDistr(0x0),
121 fCutRatioClsOverCrossRowsTPC(0),
122 fCutRatioSignalNOverCrossRowsTPC(0),
123 fCutMinCrossedRowsTPCPtDep(""),
124 f1CutMinNCrossedRowsTPCPtDep(0x0),
125 fUseCutGeoNcrNcl(kFALSE),
126 fDeadZoneWidth(3.),
127 fCutGeoNcrNclLength(130.),
128 fCutGeoNcrNclGeom1Pt(1.5),
129 fCutGeoNcrNclFractionNcr(0.85),
130 fCutGeoNcrNclFractionNcl(0.7),
131 fUseV0ANDSelectionOffline(kFALSE),
132 fUseTPCtrackCutsOnThisDaughter(kTRUE),
133 fApplyZcutOnSPDvtx(kFALSE)
134 {
135  //
136  // Default Constructor
137  //
138  fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
139 }
140 //--------------------------------------------------------------------------
142  AliAnalysisCuts(source),
143  fMinVtxType(source.fMinVtxType),
144  fMinVtxContr(source.fMinVtxContr),
146  fMaxVtxZ(source.fMaxVtxZ),
151  fTriggerMask(source.fTriggerMask),
153  fTriggerClass(),
154  fTrackCuts(0),
155  fnPtBins(source.fnPtBins),
157  fPtBinLimits(0),
158  fnVars(source.fnVars),
159  fVarNames(0),
160  fnVarsForOpt(source.fnVarsForOpt),
161  fVarsForOpt(0),
162  fGlobalIndex(source.fGlobalIndex),
163  fCutsRD(0),
164  fIsUpperCut(0),
165  fUsePID(source.fUsePID),
166  fUseAOD049(source.fUseAOD049),
167  fPidHF(0),
171  fUseMCVertex(source.fUseMCVertex),
173  fOptPileup(source.fOptPileup),
175  fMinDzPileup(source.fMinDzPileup),
181  fFixRefs(source.fFixRefs),
184  fMinPtCand(source.fMinPtCand),
185  fMaxPtCand(source.fMaxPtCand),
195  fKinkReject(source.fKinkReject),
198  fHistCentrDistr(0x0),
212 {
213  //
214  // Copy constructor
215  //
216  cout<<"Copy constructor"<<endl;
217  fTriggerClass[0] = source.fTriggerClass[0];
218  fTriggerClass[1] = source.fTriggerClass[1];
219  if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
220  if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
221  if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
222  if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
223  if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
224  if(source.fPidHF) SetPidHF(source.fPidHF);
225  if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
228  PrintAll();
229 
230 }
231 //--------------------------------------------------------------------------
233 {
234  //
235  // assignment operator
236  //
237  if(&source == this) return *this;
238 
239  AliAnalysisCuts::operator=(source);
240 
241  fMinVtxType=source.fMinVtxType;
242  fMinVtxContr=source.fMinVtxContr;
244  fMaxVtxZ=source.fMaxVtxZ;
249  fTriggerMask=source.fTriggerMask;
251  fTriggerClass[0]=source.fTriggerClass[0];
252  fTriggerClass[1]=source.fTriggerClass[1];
253  fnPtBins=source.fnPtBins;
255  fnVars=source.fnVars;
256  fGlobalIndex=source.fGlobalIndex;
257  fnVarsForOpt=source.fnVarsForOpt;
258  fUsePID=source.fUsePID;
259  fUseAOD049=source.fUseAOD049;
260  if(fPidHF) delete fPidHF;
261  fPidHF=new AliAODPidHF(*(source.GetPidHF()));
265  fUseMCVertex=source.fUseMCVertex;
267  fOptPileup=source.fOptPileup;
269  fMinDzPileup=source.fMinDzPileup;
275  fFixRefs=source.fFixRefs;
278  fMinPtCand=source.fMinPtCand;
279  fMaxPtCand=source.fMaxPtCand;
289  fKinkReject=source.fKinkReject;
291  if(fHistCentrDistr) delete fHistCentrDistr;
293  if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
294 
295  if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
296  if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
297  if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
298  if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
299  if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
300 
306 
315 
316  PrintAll();
317 
318  return *this;
319 }
320 //--------------------------------------------------------------------------
322  //
323  // Default Destructor
324  //
325  if(fTrackCuts) { delete fTrackCuts; fTrackCuts=0; }
326  if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
327  if(fVarNames) {delete [] fVarNames; fVarNames=0;}
328  if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
329  if(fCutsRD) {
330  delete [] fCutsRD;
331  fCutsRD=0;
332  }
333  if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
334  if(fPidHF){
335  delete fPidHF;
336  fPidHF=0;
337  }
339 
343  }
344 
345 }
346 //---------------------------------------------------------------------------
348  //
349  // Centrality selection
350  //
351  if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
352  AliWarning("Centrality estimator not valid");
353  return 3;
354  }else{
355  Float_t centvalue=GetCentrality((AliAODEvent*)event);
356  if (centvalue<-998.){//-999 if no centralityP
357  return 3;
359  return 3;
360  }else{
361  if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
362  return 2;
363  }
364  // selection to flatten the centrality distribution
365  if(fHistCentrDistr){
366  if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
367  }
368  }
369  }
370  return 0;
371 }
372 
373 
374 //-------------------------------------------------
375 void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
376  // set the histo for centrality flattening
377  // the centrality is flatten in the range minCentr,maxCentr
378  // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
379  // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
380  // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
381  // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
382 
383  if(maxCentr<minCentr){
384  AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
385  }
386 
388  fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
389  fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
390  Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
391  Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
392  fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
393  Double_t ref=0.,bincont=0.,binrefwidth=1.;
394  Int_t binref=0;
395  if(TMath::Abs(centrRef)<0.0001){
396  binref=fHistCentrDistr->GetMinimumBin();
397  binrefwidth=fHistCentrDistr->GetBinWidth(binref);
398  ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
399  }
400  else if(centrRef>0.){
401  binref=h->FindBin(centrRef);
402  if(binref<1||binref>h->GetNbinsX()){
403  AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
404  }
405  binrefwidth=fHistCentrDistr->GetBinWidth(binref);
406  ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
407  }
408  else{
409  if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
410  binref=fHistCentrDistr->GetMaximumBin();
411  binrefwidth=fHistCentrDistr->GetBinWidth(binref);
412  ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
413  }
414 
415  for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
416  if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
417  bincont=h->GetBinContent(j);
418  fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
419  fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
420  }
421  else{
422  h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
423  }
424  }
425 
426  fHistCentrDistr->SetBinContent(0,switchTRand);
427  return;
428 
429 }
430 
431 //-------------------------------------------------
433  //
434  // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
435  // Can be faster if it was required that fHistCentrDistr covers
436  // exactly the desired centrality range (e.g. part of the lines below should be done during the
437  // setting of the histo) and TH1::SetMinimum called
438  //
439 
440  if(!fHistCentrDistr) return kTRUE;
441  // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
442  // if(maxbin>fHistCentrDistr->GetNbinsX()){
443  // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
444  // }
445 
446  Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
447  Double_t bincont=fHistCentrDistr->GetBinContent(bin);
448  Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
449 
450  if(fHistCentrDistr->GetBinContent(0)<-0.9999){
451  if(gRandom->Uniform(1.)<bincont)return kTRUE;
452  return kFALSE;
453  }
454 
455  if(centDigits*100.<bincont)return kTRUE;
456  return kFALSE;
457 
458 }
459 //---------------------------------------------------------------------------
460 void AliRDHFCuts::SetupPID(AliVEvent *event) {
461  // Set the PID response object in the AliAODPidHF
462  // in case of old PID sets the TPC dE/dx BB parameterization
463 
464  if(fPidHF){
465  if(fPidHF->GetPidResponse()==0x0){
466  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
467  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
468  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
469  fPidHF->SetPidResponse(pidResp);
470  }
472  if(fPidHF->GetOldPid()) {
473 
474  Bool_t isMC=kFALSE;
475  TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
476  if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
477 
478  // pp, from LHC10d onwards
479  if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
480  event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
481  // pp, 2011 low energy run
482  if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
483  fPidHF->SetppLowEn2011(kTRUE);
484  fPidHF->SetOnePad(kFALSE);
485  }
486  // PbPb LHC10h
487  if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
488  // MC
489  if(isMC) fPidHF->SetMC(kTRUE);
490  if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
491  fPidHF->SetMClowenpp2011(kTRUE);
493  }else{
494  // check that AliPIDResponse object was properly set in case of using OADB
495  if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
496  }
497  }
498 }
499 //---------------------------------------------------------------------------
501  //
502  // Event selection
503  //
504  //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
505 
506 
507 
508  fWhyRejection=0;
510  Bool_t accept=kTRUE;
511 
512 
513  // check if it's MC
514  Bool_t isMC=kFALSE;
515  TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
516  if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
517 
518 
519  SetupPID(event);
520 
521  // trigger class
522  TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
523  // don't do for MC and for PbPb 2010 data
524  if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
525  if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
526  (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
527  fWhyRejection=5;
529  accept=kFALSE;
530  }
531  }
532 
533  // TEMPORARY FIX FOR GetEvent
534  Int_t nTracks=((AliAODEvent*)event)->GetNumberOfTracks();
535  for(Int_t itr=0; itr<nTracks; itr++){
536  AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
537  tr->SetAODEvent((AliAODEvent*)event);
538  }
539 
540  // TEMPORARY FIX FOR REFERENCES
541  // Fix references to daughter tracks
542  // if(fFixRefs) {
543  // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
544  // fixer->FixReferences((AliAODEvent*)event);
545  // delete fixer;
546  // }
547  //
548 
549 
550  // physics selection requirements
552  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
553  if(!isSelected) {
554  if(accept) fWhyRejection=7;
556  accept=kFALSE;
557  }else{
558  if(fUseOnlyOneTrigger){
559  if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
560  if(accept) fWhyRejection=7;
562  accept=kFALSE;
563  }
564  }
566  AliAODVZERO* v0data=(AliAODVZERO*)((AliAODEvent*)event)->GetVZEROData();
567  Int_t tv0a=v0data->GetV0ADecision();
568  Int_t tv0c=v0data->GetV0CDecision();
569  if(!(tv0a==1 && tv0c==1)){
570  if(accept) fWhyRejection=7;
572  accept=kFALSE;
573  }
574  }
575  }
576  }
577 
578  // centrality selection
579  if (fUseCentrality!=kCentOff) {
580  Int_t rejection=IsEventSelectedInCentrality(event);
581  Bool_t okCent=kFALSE;
582  if(rejection==0) okCent=kTRUE;
583  if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
584  if(!okCent){
585  if(accept) fWhyRejection=rejection;
588  accept=kFALSE;
589  }
590 
591  }
592 
593  // PbPb2011 outliers in tracklets vs. VZERO and centTRK vs. centV0
594  if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
596  Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
597  Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
598  Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
599  if(ntracklets<1000. && v0cent<cutval){
600  if(accept) fWhyRejection=2;
602  accept=kFALSE;
603  }
604  }
605  if(fMaxDiffTRKV0Centr>0.){
606  Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
607  Double_t trkcent=GetCentrality((AliAODEvent*)event,kCentTRK);
608  if(TMath::Abs(trkcent-v0cent)>fMaxDiffTRKV0Centr){
609  if(accept) fWhyRejection=1;
611  accept=kFALSE;
612  }
613  }
614  }
615 
616  // vertex requirements
617 
618  const AliVVertex *vertex = event->GetPrimaryVertex();
619 
620  if(!vertex){
621  accept=kFALSE;
623  }else{
624  TString title=vertex->GetTitle();
625  if(title.Contains("Z") && fMinVtxType>1){
626  accept=kFALSE;
628  }
629  else if(title.Contains("3D") && fMinVtxType>2){
630  accept=kFALSE;
632  }
633  if(vertex->GetNContributors()<fMinVtxContr){
634  accept=kFALSE;
636  }
637  }
638 
639  const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
640  if(fCutOnzVertexSPD>0){
641  if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
642  accept=kFALSE;
644  }else{
645  if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
646  // protection for events with bad reconstructed track vertex (introduced for 2011 Pb-Pb)
648  if(accept) fWhyRejection=6;
649  accept=kFALSE;
650  }
651  if(fCutOnzVertexSPD>=2 && vertex){
652  Double_t dz = vSPD->GetZ()-vertex->GetZ();
653  // cut on absolute distance between track and SPD vertex (introduced for 2011 Pb-Pb)
654  if(TMath::Abs(dz)>0.5) {
656  if(accept) fWhyRejection=0;
657  accept=kFALSE;
658  }
659  if(accept && fCutOnzVertexSPD==3){
660  // cut on nsigma distance between track and SPD vertex (for 2015 Pb-Pb)
661  double covTrc[6],covSPD[6];
662  vertex->GetCovarianceMatrix(covTrc);
663  vSPD->GetCovarianceMatrix(covSPD);
664  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
665  double errTrc = TMath::Sqrt(covTrc[5]);
666  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
667  if (TMath::Abs(dz)>0.2 || nsigTot>10 || nsigTrc>20){
668  // reject, bad reconstructed track vertex
670  if(accept) fWhyRejection=0;
671  accept=kFALSE;
672  }
673  }
674  }
675  }
676  }
677 
678  Double_t zvert = -999;
679  if(!fApplyZcutOnSPDvtx && vertex) zvert=vertex->GetZ();
680  else if(fApplyZcutOnSPDvtx) {
681  if(!vSPD || (vSPD && vSPD->GetNContributors()<1)){
682  accept=kFALSE;
684  }
685  else zvert = vSPD->GetZ();
686  }
687 
688  if(TMath::Abs(zvert)>fMaxVtxZ) {
690  if(accept) fWhyRejection=6;
691  accept=kFALSE;
692  }
693 
694  // pile-up rejection
696  Bool_t isPileup=kFALSE;
698  isPileup=event->IsPileupFromSPDInMultBins();
699  }else{
702  isPileup=event->IsPileupFromSPD(cutc,cutz,3.,2.,10.);
703  }
704  if(isPileup){
705  if(accept) fWhyRejection=1;
707  accept=kFALSE;
708  }
709  }
711  AliAnalysisUtils utils;
712  utils.SetMinPlpContribMV(fMinContrPileup); // min. multiplicity of the pile-up vertex to consider
713  utils.SetMaxPlpChi2MV(fMaxVtxChi2PileupMV); // max chi2 per contributor of the pile-up vertex to consider.
714  utils.SetMinWDistMV(fMinWDzPileupMV); // minimum weighted distance in Z between 2 vertices (i.e. (zv1-zv2)/sqrt(sigZv1^2+sigZv2^2) )
715  utils.SetCheckPlpFromDifferentBCMV(fRejectPlpFromDiffBCMV); // vertex with |BCID|>2 will trigger pile-up (based on TOF)
716  Bool_t isPUMV = utils.IsPileUpMV(event);
717  if(isPUMV) {
718  if(accept) fWhyRejection=1;
720  accept=kFALSE;
721  }
722  }
723 
724 
725  // Correcting PP2012 flag to remoce tracks crossing SPD misaligned staves for periods 12def
726  if(fApplySPDMisalignedPP2012 && !(event->GetRunNumber()>=195681 && event->GetRunNumber()<=197388)) fApplySPDMisalignedPP2012=false;
727 
728  return accept;
729 }
730 //---------------------------------------------------------------------------
732  //
733  // Daughter tracks selection
734  //
735  if(!fTrackCuts) return kTRUE;
736 
737  Int_t ndaughters = d->GetNDaughters();
738  AliAODVertex *vAOD = d->GetPrimaryVtx();
739  Double_t pos[3],cov[6];
740  vAOD->GetXYZ(pos);
741  vAOD->GetCovarianceMatrix(cov);
742  const AliESDVertex vESD(pos,cov,100.,100);
743 
744  Bool_t retval=kTRUE;
745 
746  for(Int_t idg=0; idg<ndaughters; idg++) {
747  AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
748  if(!dgTrack) {retval = kFALSE; continue;}
749  //printf("charge %d\n",dgTrack->Charge());
750  if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
751 
753  { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
754 
755  if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts,aod)) retval = kFALSE;
756  }
757 
758  return retval;
759 }
760 //---------------------------------------------------------------------------
762  //
763  // Check if AOD and deltaAOD files are composed of the same events:
764  // mismatches were observed in the merged AODs of LHC15o
765  //
766  // When AOD+deltaAOD are produced from ESD, mismatches can be found looking at:
767  // - the AOD trees in AliAOD.root and AliAOD.VertexingHF.root have different number of entries
768  // - the titles of the TProcessID objects do not match
769  // When deltaAOD are produced from AOD, mismatches can be found looking at:
770  // - the AOD trees in AliAOD.root and AliAOD.VertexingHF.root have different number of entries
771  //
772  // Return values:
773  // -1: AOD and deltaAOD trees have different number of entries
774  // 0: AOD and deltaAOD trees have same number of entries + the titles of the TProcessID objects do not match
775  // 1: AOD and deltaAOD trees have same number of entries + the titles of the TProcessID objects match
776  Bool_t okTProcessNames = kTRUE;
777  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
778  TTree *treeAOD = aodHandler->GetTree();
779  TTree *treeDeltaAOD = treeAOD->GetFriend("aodTree");
780  if(!treeDeltaAOD || !treeAOD) return -1;
781  if(treeDeltaAOD && treeAOD){
782  if(treeAOD->GetEntries()!=treeDeltaAOD->GetEntries()){
783  printf("AliRDHFCuts::CheckMatchingAODdeltaAODevents: Difference in number of entries in main and friend tree, skipping event\n");
784  return -1;
785  }
786  TFile *mfile = treeAOD->GetCurrentFile();
787  TFile *dfile = treeDeltaAOD->GetCurrentFile();
788  TList* lm=mfile->GetListOfKeys();
789  TList* ld=dfile->GetListOfKeys();
790  Int_t nentm=lm->GetEntries();
791  for(Int_t jm=0; jm<nentm; jm++){
792  TKey* o=(TKey*)lm->At(jm);
793  TString clnam=o->GetClassName();
794  if(clnam=="TProcessID"){
795  TString pname=o->GetName();
796  TString ptit=o->GetTitle();
797  if(pname.Contains("ProcessID")){
798  TObject* od=(TObject*)ld->FindObject(pname.Data());
799  if(od){
800  TString ptit2=od->GetTitle();
801  if(ptit2!=ptit){
802  printf("AliRDHFCuts::CheckMatchingAODdeltaAODevents: mismatch in %s: AOD: %s -- deltaAOD: %s\n",pname.Data(),ptit.Data(),ptit2.Data());
803  okTProcessNames = kFALSE;
804  }
805  }
806  }
807  }
808  }
809  }
810 
811  if (okTProcessNames) return 1;
812  else return 0;
813 }
814 //---------------------------------------------------------------------------
816  //
817  // Check the correctness of the string syntax
818  //
819  Bool_t retval=kTRUE;
820 
821  if(!rows.Contains("pt")) {
822  if(print) AliError("string must contain \"pt\"");
823  retval= kFALSE;
824  }
825  return retval;
826 }
827 //---------------------------------------------------------------------------
829  //
830  //Create the TFormula from TString for TPC crossed rows pT dependent cut
831  //
832 
833 
834  // setting data member that describes the TPC crossed rows pT dependent cut
836 
837  // creating TFormula from TString
840  // resetting TFormula
842  }
843  if(!CheckPtDepCrossedRows(rows,kTRUE))return;
844 
845  TString tmp(rows);
846  tmp.ReplaceAll("pt","x");
847  f1CutMinNCrossedRowsTPCPtDep = new TFormula("f1CutMinNCrossedRowsTPCPtDep",tmp.Data());
848 
849 
850 }
851 //---------------------------------------------------------------------------
852 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts, const AliAODEvent* aod) const{
853  //
854  // Convert to ESDtrack, relate to vertex and check cuts
855  //
856  if(!cuts) return kTRUE;
857 
858  if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
859 
860 
861  // convert to ESD track here
862  AliESDtrack esdTrack(track);
863  // set the TPC cluster info
864  esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
865  esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
866  esdTrack.SetTPCPointsF(track->GetTPCNclsF());
867  // needed to calculate the impact parameters
868  esdTrack.RelateToVertex(primary,0.,3.);
869 
870  //applying ESDtrackCut
871  if(!cuts->IsSelected(&esdTrack)) return kFALSE;
872 
873  //appliyng kink rejection
874  if(fKinkReject){
875  AliAODVertex *maybeKink=track->GetProdVertex();
876  if(maybeKink->GetType()==AliAODVertex::kKink) return kFALSE;
877  }
878 
879  //appliyng TPC crossed rows pT dependent cut
881  Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
882  if(nCrossedRowsTPC<f1CutMinNCrossedRowsTPCPtDep->Eval(esdTrack.Pt())) return kFALSE;
883  }
884 
885  //appliyng NTPCcls/NTPCcrossedRows cut
887  Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
888  Float_t nClustersTPC = esdTrack.GetTPCNcls();
889  if(nCrossedRowsTPC!=0){
890  Float_t ratio = nClustersTPC/nCrossedRowsTPC;
891  if(ratio<fCutRatioClsOverCrossRowsTPC) return kFALSE;
892  }
893  else return kFALSE;
894  }
895 
896  //appliyng TPCsignalN/NTPCcrossedRows cut
898  Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
899  Float_t nTPCsignal = esdTrack.GetTPCsignalN();
900  if(nCrossedRowsTPC!=0){
901  Float_t ratio = nTPCsignal/nCrossedRowsTPC;
902  if(ratio<fCutRatioSignalNOverCrossRowsTPC) return kFALSE;
903  }
904  else return kFALSE;
905  }
906 
907  // geometrical cut (note uses track at vertex instead of at TPC inner wall)
909  Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
910  Float_t lengthInActiveZoneTPC=esdTrack.GetLengthInActiveZone(0,fDeadZoneWidth,220.,aod->GetMagneticField());
911  Double_t cutGeoNcrNclLength=fCutGeoNcrNclLength-TMath::Power(TMath::Abs(esdTrack.GetSigned1Pt()),fCutGeoNcrNclGeom1Pt);
912  Bool_t isOK=kTRUE;
913  if (lengthInActiveZoneTPC<cutGeoNcrNclLength) isOK=kFALSE;
914  if (nCrossedRowsTPC<fCutGeoNcrNclFractionNcr*cutGeoNcrNclLength) isOK=kFALSE;
915  if (esdTrack.GetTPCncls()<fCutGeoNcrNclFractionNcl*cutGeoNcrNclLength) isOK=kFALSE;
916  if(!isOK) return kFALSE;
917  }
918 
919 
921  // to be implemented
922  // we need either to have here the AOD Event,
923  // or to have the pileup vertex object
924  }
925 
927  Bool_t deadSPDLay1PbPb2011[20][4]={
928  {kTRUE,kTRUE,kTRUE,kTRUE},
929  {kTRUE,kTRUE,kTRUE,kTRUE},
930  {kTRUE,kTRUE,kTRUE,kTRUE},
931  {kTRUE,kTRUE,kTRUE,kTRUE},
932  {kTRUE,kTRUE,kTRUE,kTRUE},
933  {kFALSE,kFALSE,kTRUE,kTRUE},
934  {kTRUE,kTRUE,kFALSE,kFALSE},
935  {kTRUE,kTRUE,kTRUE,kTRUE},
936  {kFALSE,kFALSE,kTRUE,kTRUE},
937  {kTRUE,kTRUE,kTRUE,kTRUE},
938  {kTRUE,kTRUE,kFALSE,kFALSE},
939  {kTRUE,kTRUE,kTRUE,kTRUE},
940  {kFALSE,kFALSE,kFALSE,kFALSE},
941  {kFALSE,kFALSE,kTRUE,kTRUE},
942  {kFALSE,kFALSE,kFALSE,kFALSE},
943  {kFALSE,kFALSE,kFALSE,kFALSE},
944  {kTRUE,kTRUE,kTRUE,kTRUE},
945  {kTRUE,kTRUE,kFALSE,kFALSE},
946  {kFALSE,kFALSE,kFALSE,kFALSE},
947  {kFALSE,kFALSE,kFALSE,kFALSE}
948  };
949  Bool_t deadSPDLay2PbPb2011[40][4]={
950  {kTRUE,kTRUE,kTRUE,kTRUE},
951  {kTRUE,kTRUE,kTRUE,kTRUE},
952  {kTRUE,kTRUE,kTRUE,kTRUE},
953  {kTRUE,kTRUE,kTRUE,kTRUE},
954  {kTRUE,kTRUE,kTRUE,kTRUE},
955  {kTRUE,kTRUE,kTRUE,kTRUE},
956  {kTRUE,kTRUE,kTRUE,kTRUE},
957  {kTRUE,kTRUE,kTRUE,kTRUE},
958  {kTRUE,kTRUE,kTRUE,kTRUE},
959  {kTRUE,kTRUE,kTRUE,kTRUE},
960  {kTRUE,kTRUE,kTRUE,kTRUE},
961  {kTRUE,kTRUE,kTRUE,kTRUE},
962  {kFALSE,kFALSE,kFALSE,kFALSE},
963  {kFALSE,kFALSE,kTRUE,kTRUE},
964  {kTRUE,kTRUE,kTRUE,kTRUE},
965  {kTRUE,kTRUE,kTRUE,kTRUE},
966  {kTRUE,kTRUE,kFALSE,kFALSE},
967  {kTRUE,kTRUE,kTRUE,kTRUE},
968  {kTRUE,kTRUE,kTRUE,kTRUE},
969  {kTRUE,kTRUE,kTRUE,kTRUE},
970  {kFALSE,kFALSE,kFALSE,kFALSE},
971  {kFALSE,kFALSE,kFALSE,kFALSE},
972  {kTRUE,kTRUE,kTRUE,kTRUE},
973  {kTRUE,kTRUE,kTRUE,kTRUE},
974  {kFALSE,kFALSE,kFALSE,kFALSE},
975  {kFALSE,kFALSE,kFALSE,kFALSE},
976  {kTRUE,kTRUE,kTRUE,kTRUE},
977  {kTRUE,kTRUE,kTRUE,kTRUE},
978  {kFALSE,kFALSE,kFALSE,kFALSE},
979  {kFALSE,kFALSE,kFALSE,kFALSE},
980  {kFALSE,kFALSE,kFALSE,kFALSE},
981  {kFALSE,kFALSE,kFALSE,kFALSE},
982  {kTRUE,kTRUE,kTRUE,kTRUE},
983  {kTRUE,kTRUE,kTRUE,kTRUE},
984  {kTRUE,kTRUE,kFALSE,kFALSE},
985  {kTRUE,kTRUE,kTRUE,kTRUE},
986  {kFALSE,kFALSE,kFALSE,kFALSE},
987  {kFALSE,kFALSE,kFALSE,kFALSE},
988  {kFALSE,kFALSE,kFALSE,kFALSE},
989  {kFALSE,kFALSE,kFALSE,kFALSE}
990  };
991  Double_t xyz1[3],xyz2[3];
992  esdTrack.GetXYZAt(3.9,0.,xyz1);
993  esdTrack.GetXYZAt(7.6,0.,xyz2);
994  Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
995  if(phi1<0) phi1+=2*TMath::Pi();
996  Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
997  Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
998  if(phi2<0) phi2+=2*TMath::Pi();
999  Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
1000  Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
1001  Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
1002  Bool_t lay1ok=kFALSE;
1003  if(mod1>=0 && mod1<4 && lad1<20){
1004  lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
1005  }
1006  Bool_t lay2ok=kFALSE;
1007  if(mod2>=0 && mod2<4 && lad2<40){
1008  lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
1009  }
1010  if(!lay1ok && !lay2ok) return kFALSE;
1011  }
1012 
1014  // Cut tracks crossing the SPD at 5.6<phi<2pi
1015  Double_t xyz1[3],xyz2[3];
1016  esdTrack.GetXYZAt(3.9,0.,xyz1);
1017  esdTrack.GetXYZAt(7.6,0.,xyz2);
1018  Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
1019  if(phi1<0) phi1+=2*TMath::Pi();
1020  Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
1021  if(phi2<0) phi2+=2*TMath::Pi();
1022  Bool_t lay1ok=kTRUE;
1023  if(phi1>5.6 && phi1<2.*TMath::Pi()) lay1ok=kFALSE;
1024  Bool_t lay2ok=kTRUE;
1025  if(phi2>5.6 && phi2<2.*TMath::Pi()) lay2ok=kFALSE;
1026  if(!lay1ok || !lay2ok) return kFALSE;
1027  }
1028 
1029  return kTRUE;
1030 }
1031 //---------------------------------------------------------------------------
1032 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
1033  // Set the pt bins
1034 
1035  if(fPtBinLimits) {
1036  delete [] fPtBinLimits;
1037  fPtBinLimits = NULL;
1038  printf("Changing the pt bins\n");
1039  }
1040 
1041  if(nPtBinLimits != fnPtBins+1){
1042  cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
1043  SetNPtBins(nPtBinLimits-1);
1044  }
1045 
1046  fnPtBinLimits = nPtBinLimits;
1047  SetGlobalIndex();
1048  //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
1050  for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
1051 
1052  return;
1053 }
1054 //---------------------------------------------------------------------------
1055 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
1056  // Set the variable names
1057 
1058  if(fVarNames) {
1059  delete [] fVarNames;
1060  fVarNames = NULL;
1061  //printf("Changing the variable names\n");
1062  }
1063  if(nVars!=fnVars){
1064  printf("Wrong number of variables: it has to be %d\n",fnVars);
1065  return;
1066  }
1067  //fnVars=nVars;
1068  fVarNames = new TString[nVars];
1069  fIsUpperCut = new Bool_t[nVars];
1070  for(Int_t iv=0; iv<nVars; iv++) {
1071  fVarNames[iv] = varNames[iv];
1072  fIsUpperCut[iv] = isUpperCut[iv];
1073  }
1074 
1075  return;
1076 }
1077 //---------------------------------------------------------------------------
1079  // Set the variables to be used for cuts optimization
1080 
1081  if(fVarsForOpt) {
1082  delete [] fVarsForOpt;
1083  fVarsForOpt = NULL;
1084  //printf("Changing the variables for cut optimization\n");
1085  }
1086 
1087  if(nVars==0){
1088  printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
1089  return;
1090  }
1091 
1092  fnVarsForOpt = 0;
1093  fVarsForOpt = new Bool_t[fnVars];
1094  for(Int_t iv=0; iv<fnVars; iv++) {
1095  fVarsForOpt[iv]=forOpt[iv];
1096  if(fVarsForOpt[iv]) fnVarsForOpt++;
1097  }
1098 
1099  return;
1100 }
1101 
1102 //---------------------------------------------------------------------------
1104  //
1105  // set centrality estimator
1106  //
1107  fUseCentrality=flag;
1108  if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
1109 
1110  return;
1111 }
1112 
1113 
1114 //---------------------------------------------------------------------------
1116  //
1117  // store the cuts
1118  //
1119  if(nVars!=fnVars) {
1120  printf("Wrong number of variables: it has to be %d\n",fnVars);
1121  AliFatal("exiting");
1122  }
1123  if(nPtBins!=fnPtBins) {
1124  printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
1125  AliFatal("exiting");
1126  }
1127 
1128  if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
1129 
1130 
1131  for(Int_t iv=0; iv<fnVars; iv++) {
1132 
1133  for(Int_t ib=0; ib<fnPtBins; ib++) {
1134 
1135  //check
1136  if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
1137  cout<<"Overflow, exit..."<<endl;
1138  return;
1139  }
1140 
1141  fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
1142 
1143  }
1144  }
1145  return;
1146 }
1147 //---------------------------------------------------------------------------
1148 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
1149  //
1150  // store the cuts
1151  //
1152  if(glIndex != fGlobalIndex){
1153  cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
1154  AliFatal("exiting");
1155  }
1156  if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
1157 
1158  for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
1159  fCutsRD[iGl] = cutsRDGlob[iGl];
1160  }
1161  return;
1162 }
1163 //---------------------------------------------------------------------------
1165  //
1166  // print all cuts values
1167  //
1168 
1169  printf("Minimum vtx type %d\n",fMinVtxType);
1170  printf("Minimum vtx contr %d\n",fMinVtxContr);
1171  printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
1172  printf("Min SPD mult %d\n",fMinSPDMultiplicity);
1173  printf("Use PID %d OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
1174  printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
1175  printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
1176  printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
1177  if(fOptPileup==1) printf(" -- Reject pileup event");
1178  if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
1179  if(fUseCentrality>0) {
1180  TString estimator="";
1181  if(fUseCentrality==kCentV0M) estimator = "V0";
1182  if(fUseCentrality==kCentTRK) estimator = "Tracks";
1183  if(fUseCentrality==kCentTKL) estimator = "Tracklets";
1184  if(fUseCentrality==kCentCL1) estimator = "SPD clusters outer";
1185  if(fUseCentrality==kCentZNA) estimator = "ZNA";
1186  if(fUseCentrality==kCentZPA) estimator = "ZPA";
1187  if(fUseCentrality==kCentV0A) estimator = "V0A";
1188  if(fUseCentrality==kCentCL0) estimator = "SPD clusters inner";
1189  printf("Centrality class considered: %.1f-%.1f, estimated with %s\n",fMinCentrality,fMaxCentrality,estimator.Data());
1190  }
1191  if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
1192 
1193  if(fCutRatioClsOverCrossRowsTPC) printf("N TPC Clusters > %f N TPC Crossed Rows\n", fCutRatioClsOverCrossRowsTPC);
1194  if(fCutRatioSignalNOverCrossRowsTPC) printf("N TPC Points for dE/dx > %f N TPC Crossed Rows\n", fCutRatioSignalNOverCrossRowsTPC);
1195  if(f1CutMinNCrossedRowsTPCPtDep) printf("N TPC Crossed Rows pT-dependent cut: %s\n", fCutMinCrossedRowsTPCPtDep.Data());
1196 
1197  if(fVarNames){
1198  cout<<"Array of variables"<<endl;
1199  for(Int_t iv=0;iv<fnVars;iv++){
1200  cout<<fVarNames[iv]<<"\t";
1201  }
1202  cout<<endl;
1203  }
1204  if(fVarsForOpt){
1205  cout<<"Array of optimization"<<endl;
1206  for(Int_t iv=0;iv<fnVars;iv++){
1207  cout<<fVarsForOpt[iv]<<"\t";
1208  }
1209  cout<<endl;
1210  }
1211  if(fIsUpperCut){
1212  cout<<"Array of upper/lower cut"<<endl;
1213  for(Int_t iv=0;iv<fnVars;iv++){
1214  cout<<fIsUpperCut[iv]<<"\t";
1215  }
1216  cout<<endl;
1217  }
1218  if(fPtBinLimits){
1219  cout<<"Array of ptbin limits"<<endl;
1220  for(Int_t ib=0;ib<fnPtBinLimits;ib++){
1221  cout<<fPtBinLimits[ib]<<"\t";
1222  }
1223  cout<<endl;
1224  }
1225  if(fCutsRD){
1226  cout<<"Matrix of cuts"<<endl;
1227  for(Int_t iv=0;iv<fnVars;iv++){
1228  for(Int_t ib=0;ib<fnPtBins;ib++){
1229  cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
1230  }
1231  cout<<endl;
1232  }
1233  cout<<endl;
1234  }
1235  if(fPidHF) fPidHF->PrintAll();
1236  return;
1237 }
1238 
1239 //--------------------------------------------------------------------------
1241  // print the trigger selection
1242 
1243  printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
1244 
1245  cout<<" Trigger selection pattern: ";
1246  if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
1247  if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
1248  if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
1249  if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
1250  if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
1251  if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
1252  if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
1253  if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
1254  if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
1255  if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
1256  cout << endl<< endl;
1257 
1258 }
1259 
1260 //---------------------------------------------------------------------------
1261 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
1262  //
1263  // get the cuts
1264  //
1265 
1266  //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
1267 
1268 
1269  Int_t iv,ib;
1270  if(!cutsRD) {
1271  //cout<<"Initialization..."<<endl;
1272  cutsRD=new Float_t*[fnVars];
1273  for(iv=0; iv<fnVars; iv++) {
1274  cutsRD[iv] = new Float_t[fnPtBins];
1275  }
1276  }
1277 
1278  for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
1279  GetVarPtIndex(iGlobal,iv,ib);
1280  cutsRD[iv][ib] = fCutsRD[iGlobal];
1281  }
1282 
1283  return;
1284 }
1285 
1286 //---------------------------------------------------------------------------
1288  //
1289  // give the global index from variable and pt bin
1290  //
1291  return iPtBin*fnVars+iVar;
1292 }
1293 
1294 //---------------------------------------------------------------------------
1295 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
1296  //
1297  //give the index of the variable and of the pt bin from the global index
1298  //
1299  iPtBin=(Int_t)iGlob/fnVars;
1300  iVar=iGlob%fnVars;
1301 
1302  return;
1303 }
1304 
1305 //---------------------------------------------------------------------------
1307  //
1308  //give the pt bin where the pt lies.
1309  //
1310  Int_t ptbin=-1;
1311  if(pt<fPtBinLimits[0])return ptbin;
1312  for (Int_t i=0;i<fnPtBins;i++){
1313  if(pt<fPtBinLimits[i+1]) {
1314  ptbin=i;
1315  break;
1316  }
1317  }
1318  return ptbin;
1319 }
1320 //-------------------------------------------------------------------
1322  //
1323  // Give the value of cut set for the variable iVar and the pt bin iPtBin
1324  //
1325  if(!fCutsRD){
1326  cout<<"Cuts not iniziaisez yet"<<endl;
1327  return 0;
1328  }
1329  return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1330 }
1331 
1332 //-------------------------------------------------------------------
1334 
1335  if(aodEvent->GetRunNumber()<244824)return GetCentralityOldFramework(aodEvent,estimator);
1336  Double_t cent=-999;
1337 
1338  if(estimator==kCentOff) return -999;
1339 
1340  AliMultSelection *multSelection = (AliMultSelection*)aodEvent->FindListObject(fMultSelectionObjectName);
1341  if(!multSelection){
1342  AliWarning("AliMultSelection could not be found in the aod event list of objects");
1343  return cent;
1344  }
1345 
1346  // Compare centrality with the centrality at AOD filtering and on-the-fly
1347  if( fMultSelectionObjectName.CompareTo("MultSelection")!=0 ){
1348  AliMultSelection *defaultmultSelection = (AliMultSelection*)aodEvent->FindListObject("MultSelection");
1349  if(!defaultmultSelection){
1350  AliWarning("AliMultSelection default method could not be found in the aod event list of objects");
1351  return cent;
1352  }
1353  Float_t defaultCent = defaultmultSelection->GetMultiplicityPercentile("V0M");
1354  Float_t newCent = multSelection->GetMultiplicityPercentile("V0M");
1355  if( defaultCent<20. && newCent>20.) fEvRejectionBits+=1<<kMismatchOldNewCentrality;
1356  else if (defaultCent>20. && newCent<20.) fEvRejectionBits+=1<<kMismatchOldNewCentrality;
1357  }
1358 
1359  if(estimator==kCentV0M){
1360  cent=multSelection->GetMultiplicityPercentile("V0M");
1361  }else if(estimator==kCentV0A){
1362  cent=multSelection->GetMultiplicityPercentile("V0A");
1363  }else if(estimator==kCentZNA){
1364  cent=multSelection->GetMultiplicityPercentile("ZNA");
1365  }else if(estimator==kCentCL1){
1366  cent=multSelection->GetMultiplicityPercentile("CL1");
1367  }else if(estimator==kCentCL0){
1368  cent=multSelection->GetMultiplicityPercentile("CL0");
1369  }else {
1370  AliWarning(Form("CENTRALITY ESTIMATE WITH ESTIMATOR %d NOT YET IMPLEMENTED FOR NEW FRAMEWORK",(Int_t)estimator));
1371  return cent;
1372  }
1373  Int_t qual = multSelection->GetEvSelCode();
1374  if(qual == 199 ) cent=-999;
1375  return cent;
1376 }
1377 //-------------------------------------------------------------------
1379  //
1380  // Get centrality percentile
1381  //
1382 
1383  TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1384  if(mcArray) {fUseAOD049=kFALSE;}
1385 
1386  AliAODHeader *header=dynamic_cast<AliAODHeader*>(aodEvent->GetHeader());
1387  if(!header) AliFatal("Not a standard AOD");
1388  AliCentrality *centrality=header->GetCentralityP();
1389  Float_t cent=-999.;
1390  Bool_t isSelRun=kFALSE;
1391  Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1392  if(!centrality) return cent;
1393  else{
1394  if (estimator==kCentV0M){
1395  cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1396  if(cent<0){
1397  Int_t quality = centrality->GetQuality();
1398  if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1399  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1400  }else{
1401  Int_t runnum=aodEvent->GetRunNumber();
1402  for(Int_t ir=0;ir<5;ir++){
1403  if(runnum==selRun[ir]){
1404  isSelRun=kTRUE;
1405  break;
1406  }
1407  }
1408  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1409  }
1410  }
1411 
1412  //temporary fix for AOD049 outliers
1413  if(fUseAOD049&&cent>=0){
1414  Float_t v0=0;
1415  AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1416  v0+=aodV0->GetMTotV0A();
1417  v0+=aodV0->GetMTotV0C();
1418  if(cent==0&&v0<19500)return -1;//filtering issue
1419  Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1420  Float_t val= 1.30552 + 0.147931 * v0;
1421  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};
1422  if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1423  }
1424  }
1425  else {
1426  if (estimator==kCentTRK) {
1427  cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1428  if(cent<0){
1429  Int_t quality = centrality->GetQuality();
1430  if(quality<=1){
1431  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1432  }else{
1433  Int_t runnum=aodEvent->GetRunNumber();
1434  for(Int_t ir=0;ir<5;ir++){
1435  if(runnum==selRun[ir]){
1436  isSelRun=kTRUE;
1437  break;
1438  }
1439  }
1440  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1441  }
1442  }
1443  }
1444  else{
1445  if (estimator==kCentTKL){
1446  cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1447  if(cent<0){
1448  Int_t quality = centrality->GetQuality();
1449  if(quality<=1){
1450  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1451  }else{
1452  Int_t runnum=aodEvent->GetRunNumber();
1453  for(Int_t ir=0;ir<5;ir++){
1454  if(runnum==selRun[ir]){
1455  isSelRun=kTRUE;
1456  break;
1457  }
1458  }
1459  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1460  }
1461  }
1462  }
1463  else{
1464  if (estimator==kCentCL1){
1465  cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1466  if(cent<0){
1467  Int_t quality = centrality->GetQuality();
1468  if(quality<=1){
1469  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1470  }else{
1471  Int_t runnum=aodEvent->GetRunNumber();
1472  for(Int_t ir=0;ir<5;ir++){
1473  if(runnum==selRun[ir]){
1474  isSelRun=kTRUE;
1475  break;
1476  }
1477  }
1478  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1479  }
1480  }
1481  }
1482  else{
1483  if (estimator==kCentZNA){
1484  cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1485  if(cent<0){
1486  Int_t quality = centrality->GetQuality();
1487  if(quality<=1){
1488  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1489  }else{
1490  Int_t runnum=aodEvent->GetRunNumber();
1491  for(Int_t ir=0;ir<5;ir++){
1492  if(runnum==selRun[ir]){
1493  isSelRun=kTRUE;
1494  break;
1495  }
1496  }
1497  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1498  }
1499  }
1500  }
1501  else{
1502  if (estimator==kCentZPA){
1503  cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1504  if(cent<0){
1505  Int_t quality = centrality->GetQuality();
1506  if(quality<=1){
1507  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1508  }else{
1509  Int_t runnum=aodEvent->GetRunNumber();
1510  for(Int_t ir=0;ir<5;ir++){
1511  if(runnum==selRun[ir]){
1512  isSelRun=kTRUE;
1513  break;
1514  }
1515  }
1516  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1517  }
1518  }
1519  }
1520  else{
1521  if (estimator==kCentV0A){
1522  cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1523  if(cent<0){
1524  Int_t quality = centrality->GetQuality();
1525  if(quality<=1){
1526  cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1527  }else{
1528  Int_t runnum=aodEvent->GetRunNumber();
1529  for(Int_t ir=0;ir<5;ir++){
1530  if(runnum==selRun[ir]){
1531  isSelRun=kTRUE;
1532  break;
1533  }
1534  }
1535  if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1536  }
1537  }
1538  }
1539  else {
1540  AliWarning("Centrality estimator not valid");
1541 
1542  }
1543  }
1544  }
1545  }
1546  }
1547  }
1548  }
1549  }
1550  return cent;
1551 }
1552 //-------------------------------------------------------------------
1554  //
1555  // Compare two cuts objects
1556  //
1557 
1558  Bool_t areEqual=kTRUE;
1559 
1560  if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1561 
1562  if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1563 
1564  if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1565 
1566  if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1567 
1568  if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1569 
1570  if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1571  if(fTrackCuts){
1572  if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1573 
1574  if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1575 
1576  if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1577 
1578  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;}
1579  }
1580 
1581  if(fCutsRD) {
1582  for(Int_t iv=0;iv<fnVars;iv++) {
1583  for(Int_t ib=0;ib<fnPtBins;ib++) {
1584  if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1585  cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1586  areEqual=kFALSE;
1587  }
1588  }
1589  }
1590  }
1591 
1592  return areEqual;
1593 }
1594 //---------------------------------------------------------------------------
1596  //
1597  // print cuts values in table format
1598  //
1599 
1600  TString ptString = "pT range";
1601  if(fVarNames && fPtBinLimits && fCutsRD){
1602  TString firstLine(Form("* %-15s",ptString.Data()));
1603  for (Int_t ivar=0; ivar<fnVars; ivar++){
1604  firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1605  if (ivar == fnVars){
1606  firstLine+="*\n";
1607  }
1608  }
1609  Printf("%s",firstLine.Data());
1610 
1611  for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1612  TString line;
1613  if (ipt==fnPtBins-1){
1614  line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1615  }
1616  else{
1617  line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1618  }
1619  for (Int_t ivar=0; ivar<fnVars; ivar++){
1620  line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1621  }
1622  Printf("%s",line.Data());
1623  }
1624 
1625  }
1626 
1627 
1628  return;
1629 }
1630 //--------------------------------------------------------------------------
1632  AliAODEvent *aod) const
1633 {
1634  //
1635  // Recalculate primary vertex without daughters
1636  //
1637 
1638  if(!aod) {
1639  AliError("Can not remove daughters from vertex without AOD event");
1640  return 0;
1641  }
1642 
1643  AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1644  if(!recvtx){
1645  AliDebug(2,"Removal of daughter tracks failed");
1646  return kFALSE;
1647  }
1648 
1649 
1650  //set recalculed primary vertex
1651  d->SetOwnPrimaryVtx(recvtx);
1652  delete recvtx;
1653 
1654  return kTRUE;
1655 }
1656 //--------------------------------------------------------------------------
1658 {
1659  //
1660  // Recalculate primary vertex without daughters
1661  //
1662 
1663  if(!aod) {
1664  AliError("Can not get MC vertex without AOD event");
1665  return kFALSE;
1666  }
1667 
1668  // load MC header
1669  AliAODMCHeader *mcHeader =
1670  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1671  if(!mcHeader) {
1672  AliError("Can not get MC vertex without AODMCHeader event");
1673  return kFALSE;
1674  }
1675  Double_t pos[3];
1676  Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1677  mcHeader->GetVertex(pos);
1678  AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1679 
1680  if(!recvtx){
1681  AliDebug(2,"Removal of daughter tracks failed");
1682  return kFALSE;
1683  }
1684 
1685  //set recalculed primary vertex
1686  d->SetOwnPrimaryVtx(recvtx);
1687 
1688  d->RecalculateImpPars(recvtx,aod);
1689 
1690  delete recvtx;
1691 
1692  return kTRUE;
1693 }
1694 //--------------------------------------------------------------------------
1696  AliAODEvent *aod,
1697  AliAODVertex *origownvtx) const
1698 {
1699  //
1700  // Clean-up own primary vertex if needed
1701  //
1702 
1704  d->UnsetOwnPrimaryVtx();
1705  if(origownvtx) {
1706  d->SetOwnPrimaryVtx(origownvtx);
1707  delete origownvtx; origownvtx=NULL;
1708  }
1709  d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1710  } else {
1711  if(origownvtx) {
1712  delete origownvtx; origownvtx=NULL;
1713  }
1714  }
1715  return;
1716 }
1717 //--------------------------------------------------------------------------
1718 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1719 {
1720  //
1721  // Checks if this candidate is matched to MC signal
1722  //
1723 
1724  if(!aod) return kFALSE;
1725 
1726  // get the MC array
1727  TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1728 
1729  if(!mcArray) return kFALSE;
1730 
1731  // try to match
1732  Int_t label = d->MatchToMC(pdg,mcArray);
1733 
1734  if(label>=0) {
1735  //printf("MATCH!\n");
1736  return kTRUE;
1737  }
1738 
1739  return kFALSE;
1740 }
1741 
1742 
1743 //--------------------------------------------------------------------------
1745  // recompute event primary vertex from AOD tracks
1746 
1747  AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1748  vertexer->SetITSMode();
1749  vertexer->SetMinClusters(3);
1750 
1751  AliAODVertex* pvtx=event->GetPrimaryVertex();
1752  if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1753  Float_t diamondcovxy[3];
1754  event->GetDiamondCovXY(diamondcovxy);
1755  Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1756  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1757  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1758  vertexer->SetVtxStart(diamond);
1759  delete diamond; diamond=NULL;
1760  }
1761 
1762  AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1763  if(!vertexESD) return kFALSE;
1764  if(vertexESD->GetNContributors()<=0) {
1765  //AliDebug(2,"vertexing failed");
1766  delete vertexESD; vertexESD=NULL;
1767  return kFALSE;
1768  }
1769  delete vertexer; vertexer=NULL;
1770 
1771  // convert to AliAODVertex
1772  Double_t pos[3],cov[6],chi2perNDF;
1773  vertexESD->GetXYZ(pos); // position
1774  vertexESD->GetCovMatrix(cov); //covariance matrix
1775  chi2perNDF = vertexESD->GetChi2toNDF();
1776  delete vertexESD; vertexESD=NULL;
1777 
1778  pvtx->SetPosition(pos[0],pos[1],pos[2]);
1779  pvtx->SetChi2perNDF(chi2perNDF);
1780  pvtx->SetCovMatrix(cov);
1781 
1782  return kTRUE;
1783 }
void SetNPtBins(Int_t nptBins)
Definition: AliRDHFCuts.h:403
Int_t fIsSelectedCuts
fix the daughter track references
Definition: AliRDHFCuts.h:455
Int_t pdg
Bool_t IsEventSelectedForCentrFlattening(Float_t centvalue)
Float_t fMinDzPileup
min. n. of tracklets in pileup vertex
Definition: AliRDHFCuts.h:444
Double_t fMaxPtCandTrackSPDFirst
flag to select the track kFirst criteria for pt < ptlimit
Definition: AliRDHFCuts.h:462
Int_t fWhyRejection
PID for heavy flavours manager.
Definition: AliRDHFCuts.h:437
double Double_t
Definition: External.C:58
Int_t fCutOnzVertexSPD
flag to apply cut on tracklets vs. centrality for 2011 data
Definition: AliRDHFCuts.h:467
void SetupPID(AliVEvent *event)
Bool_t IsSignalMC(AliAODRecoDecay *d, AliAODEvent *aod, Int_t pdg) const
const char * title
Definition: MakeQAPdf.C:27
Double_t fCutGeoNcrNclFractionNcr
3rd parameter of GeoNcrNcl cut
Definition: AliRDHFCuts.h:480
Bool_t fUseMCVertex
flag to switch on the removal of duaghters from the primary vertex computation
Definition: AliRDHFCuts.h:440
Bool_t fUseOnlyOneTrigger
trigger mask
Definition: AliRDHFCuts.h:419
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t * fIsUpperCut
Definition: AliRDHFCuts.h:433
Bool_t SetMCPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Bool_t fRemoveDaughtersFromPrimary
Definition: AliRDHFCuts.h:439
void SetHistoForCentralityFlattening(TH1F *h, Double_t minCentr, Double_t maxCentr, Double_t centrRef=0., Int_t switchTRand=0)
Bool_t GetUseCombined()
Definition: AliAODPidHF.h:165
Bool_t fUseMultDepPileupCut
min deltaz between main and pileup vertices
Definition: AliRDHFCuts.h:445
void SetUseCentrality(Int_t flag=1)
Float_t * fPtBinLimits
"number of limits", that is fnPtBins+1
Definition: AliRDHFCuts.h:426
Int_t fMinVtxContr
0: not cut; 1: SPDZ; 2: SPD3D; 3: Tracks
Definition: AliRDHFCuts.h:411
centrality
char Char_t
Definition: External.C:18
static Int_t CheckMatchingAODdeltaAODevents()
Float_t fMaxVtxRedChi2
minimum vertex contributors
Definition: AliRDHFCuts.h:412
void SetMClowenpp2011(Bool_t mc)
Definition: AliAODPidHF.h:102
const Float_t * GetCuts() const
Definition: AliRDHFCuts.h:254
void SetMinCrossedRowsTPCPtDep(const char *rows="")
Float_t fCutRatioClsOverCrossRowsTPC
histogram with reference centrality distribution for centrality distribution flattening ...
Definition: AliRDHFCuts.h:472
Double_t fDeadZoneWidth
flag for enabling/disabling geometrical cut on TPC track
Definition: AliRDHFCuts.h:477
Int_t fMinVtxType
cuts on the event
Definition: AliRDHFCuts.h:410
void SetppLowEn2011(Bool_t opt)
Definition: AliAODPidHF.h:104
Double_t fMaxRapidityCand
minimum pt of the candidate
Definition: AliRDHFCuts.h:459
Double_t fCutGeoNcrNclFractionNcl
4th parameter of GeoNcrNcl cut
Definition: AliRDHFCuts.h:481
void SetGlobalIndex()
Definition: AliRDHFCuts.h:201
AliRDHFCuts & operator=(const AliRDHFCuts &source)
Float_t GetCutValue(Int_t iVar, Int_t iPtBin) const
Int_t fUseCentrality
flag to use a multiplicity dependent pileup selection
Definition: AliRDHFCuts.h:446
Bool_t fUsePID
Definition: AliRDHFCuts.h:434
TRandom * gRandom
Bool_t fRejectPlpFromDiffBCMV
minimum weighted distance in Z between 2 vertices (multi-vertexer)
Definition: AliRDHFCuts.h:417
void PrintAll() const
void SetPidHF(AliAODPidHF *pidObj)
see enum below
Definition: AliRDHFCuts.h:215
Bool_t GetOldPid()
Definition: AliAODPidHF.h:157
TList * GetList(const TDirectory *d, const char *what)
Definition: ExtractData.C:172
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:247
void RecalculateImpPars(AliAODVertex *vtxAODNew, AliAODEvent *aod)
Bool_t fRemoveTrackletOutliers
Max. difference between TRK and V0 centrality (remove TPC pileup for PbPb 2011)
Definition: AliRDHFCuts.h:466
Bool_t fUseV0ANDSelectionOffline
5th parameter of GeoNcrNcl cut
Definition: AliRDHFCuts.h:482
AliRDHFCuts(const Char_t *name="RDHFCuts", const Char_t *title="")
Definition: AliRDHFCuts.cxx:64
const Int_t nPtBins
ULong64_t fTriggerMask
flag to reject pileup from different BC (multi-vertexer)
Definition: AliRDHFCuts.h:418
Bool_t fUseCutGeoNcrNcl
pT-dep cut in TPC minimum n crossed rows
Definition: AliRDHFCuts.h:476
int Int_t
Definition: External.C:63
Bool_t fUseTrackSelectionWithFilterBits
flag to reject kink daughters
Definition: AliRDHFCuts.h:469
void SetCuts(Int_t nVars, Int_t nPtBins, Float_t **cutsRD)
Int_t fnVarsForOpt
Definition: AliRDHFCuts.h:429
Float_t fMinWDzPileupMV
max chi2 per contributor of the pile-up vertex to consider (multi-vertexer).
Definition: AliRDHFCuts.h:416
Double_t fCutGeoNcrNclLength
1st parameter of GeoNcrNcl cut
Definition: AliRDHFCuts.h:478
float Float_t
Definition: External.C:68
virtual ~AliRDHFCuts()
void SetPbPb(Bool_t pbpb)
Definition: AliAODPidHF.h:105
TString * fVarNames
number of cut vars for candidates
Definition: AliRDHFCuts.h:428
Double_t fMaxPtCand
minimum pt of the candidate
Definition: AliRDHFCuts.h:458
AliESDtrackCuts * fTrackCuts
quality cuts on the daughter tracks
Definition: AliRDHFCuts.h:422
void MakeTable() const
Bool_t fFixRefs
name of the AliMultSelection object to be considered
Definition: AliRDHFCuts.h:454
Bool_t fKeepSignalMC
max rapidity of candidate (if !=-999 overrides IsInFiducialAcceptance)
Definition: AliRDHFCuts.h:460
Int_t fOptPileup
use Physics selection criteria
Definition: AliRDHFCuts.h:442
AliESDtrackCuts * GetTrackCuts() const
Definition: AliRDHFCuts.h:262
void PrintTrigger() const
Float_t fMaxVtxZ
maximum chi2/ndf
Definition: AliRDHFCuts.h:413
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:258
Int_t fIsSelectedPID
outcome of cuts selection
Definition: AliRDHFCuts.h:456
Float_t fMinCentrality
Definition: AliRDHFCuts.h:451
AliAODVertex * RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod)
void SetUpCombinedPID()
AliPIDResponse * GetPidResponse() const
Definition: AliAODPidHF.h:160
UInt_t fEvRejectionBits
used to code the step at which candidate was rejected
Definition: AliRDHFCuts.h:438
Double_t fMaxDiffTRKV0Centr
flag to apply cut on tracks crossing SPD misaligned modules for PP2012 data
Definition: AliRDHFCuts.h:465
Float_t GetCentralityOldFramework(AliAODEvent *aodEvent, AliRDHFCuts::ECentrality estimator)
Int_t fnPtBinLimits
number of pt bins for cuts
Definition: AliRDHFCuts.h:425
Float_t * fCutsRD
fnVars*fnPtBins
Definition: AliRDHFCuts.h:432
Bool_t CheckPtDepCrossedRows(TString rows, Bool_t print=kFALSE) const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Bool_t fUseCentrFlatteningInMC
flag to enable/disable the check on filter bits
Definition: AliRDHFCuts.h:470
void SetVarsForOpt(Int_t nVars, Bool_t *forOpt)
void SetVarNames(Int_t nVars, TString *varNames, Bool_t *isUpperCut)
void GetVarPtIndex(Int_t iGlob, Int_t &iVar, Int_t &iPtBin) const
void SetBetheBloch()
Bool_t * fVarsForOpt
number of cut vars to be optimized for candidates
Definition: AliRDHFCuts.h:430
Bool_t isMC
Bool_t CompareCuts(const AliRDHFCuts *obj) const
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
Bool_t fApplySPDMisalignedPP2012
flag to apply SPD dead module map of PbPb2011
Definition: AliRDHFCuts.h:464
Bool_t IsEventSelected(AliVEvent *event)
Int_t fMinSPDMultiplicity
maximum |z| of primary vertex
Definition: AliRDHFCuts.h:414
void SetMC(Bool_t mc)
Definition: AliAODPidHF.h:101
virtual void PrintAll() const
Bool_t fUseAOD049
enable PID usage (off by default)
Definition: AliRDHFCuts.h:435
Bool_t fApplySPDDeadPbPb2011
maximum pt of the candidate for which to check if the daughters fulfill kFirst criteria ...
Definition: AliRDHFCuts.h:463
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
Float_t fMaxVtxChi2PileupMV
SPD multiplicity.
Definition: AliRDHFCuts.h:415
Bool_t fIsCandTrackSPDFirst
IsSelected returns always kTRUE for MC signal.
Definition: AliRDHFCuts.h:461
TH1F * fHistCentrDistr
flag for enabling/diabling centrality flattening in MC
Definition: AliRDHFCuts.h:471
void SetPtBins(Int_t nPtBinLimits, Float_t *ptBinLimits)
Bool_t fKinkReject
cut on zSPD vertex to remove outliers in centrality vs. tracklets (0=no cut, 1= cut at 12 cm...
Definition: AliRDHFCuts.h:468
TString fMultSelectionObjectName
maximum centrality for selected events
Definition: AliRDHFCuts.h:453
Int_t fMinContrPileup
option for pielup selection
Definition: AliRDHFCuts.h:443
Bool_t fUseTPCtrackCutsOnThisDaughter
flag to apply V0AND selection offline
Definition: AliRDHFCuts.h:483
Bool_t fUsePhysicsSelection
use MC primary vertex
Definition: AliRDHFCuts.h:441
TFormula * f1CutMinNCrossedRowsTPCPtDep
pT-dep cut in TPC minimum n crossed rows
Definition: AliRDHFCuts.h:475
AliAODVertex * GetPrimaryVtx() const
void AddTrackCuts(const AliESDtrackCuts *cuts)
Definition: AliRDHFCuts.h:207
TString fTriggerClass[2]
flag to select one trigger only
Definition: AliRDHFCuts.h:420
bool Bool_t
Definition: External.C:53
Int_t fnPtBins
cuts on the candidate
Definition: AliRDHFCuts.h:424
Double_t fCutGeoNcrNclGeom1Pt
2nd parameter of GeoNcrNcl cut
Definition: AliRDHFCuts.h:479
Bool_t fApplyZcutOnSPDvtx
flag to apply TPC track quality cuts on specific D-meson daughter (used for different strategies for ...
Definition: AliRDHFCuts.h:484
AliAODPidHF * fPidHF
enable AOD049 centrality cleanup
Definition: AliRDHFCuts.h:436
Float_t fCutRatioSignalNOverCrossRowsTPC
min. value ratio NTPCClusters/NTPCCrossedRows, cut if !=0
Definition: AliRDHFCuts.h:473
Int_t fGlobalIndex
Definition: AliRDHFCuts.h:431
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
Int_t PtBin(Double_t pt) const
Int_t GetGlobalIndex(Int_t iVar, Int_t iPtBin) const
TString fCutMinCrossedRowsTPCPtDep
min. value ratio TPCPointsUsedForPID/NTPCCrossedRows, cut if !=0
Definition: AliRDHFCuts.h:474
Double_t fMinPtCand
outcome of PID selection
Definition: AliRDHFCuts.h:457
Bool_t IsDaughterSelected(AliAODTrack *track, const AliESDVertex *primary, AliESDtrackCuts *cuts, const AliAODEvent *aod=0x0) const
void SetOnePad(Bool_t onepad)
Definition: AliAODPidHF.h:103
Float_t fMaxCentrality
minimum centrality for selected events
Definition: AliRDHFCuts.h:452
void SetPidResponse(AliPIDResponse *pidResp)
Definition: AliAODPidHF.h:111
Bool_t RecomputePrimaryVertex(AliAODEvent *event) const
Int_t fnVars
Definition: AliRDHFCuts.h:427