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