1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // Base class for cuts on AOD reconstructed heavy-flavour decay
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24 #include <Riostream.h>
26 #include "AliVEvent.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliVVertex.h"
30 #include "AliESDVertex.h"
32 #include "AliAODVertex.h"
33 #include "AliESDtrack.h"
34 #include "AliAODTrack.h"
35 #include "AliESDtrackCuts.h"
36 #include "AliCentrality.h"
37 #include "AliAODRecoDecayHF.h"
38 #include "AliAnalysisVertexingHF.h"
39 #include "AliAODMCHeader.h"
40 #include "AliAODMCParticle.h"
41 #include "AliVertexerTracks.h"
42 #include "AliRDHFCuts.h"
43 #include "AliAnalysisManager.h"
44 #include "AliInputEventHandler.h"
45 #include "AliPIDResponse.h"
46 #include "AliAnalysisUtils.h"
55 //--------------------------------------------------------------------------
56 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
57 AliAnalysisCuts(name,title),
62 fMinSPDMultiplicity(0),
63 fTriggerMask(AliVEvent::kAnyINT),
64 fUseOnlyOneTrigger(kFALSE),
81 fRemoveDaughtersFromPrimary(kFALSE),
82 fRecomputePrimVertex(kFALSE),
84 fUsePhysicsSelection(kTRUE),
96 fMaxRapidityCand(-999.),
97 fKeepSignalMC(kFALSE),
98 fIsCandTrackSPDFirst(kFALSE),
99 fMaxPtCandTrackSPDFirst(0.),
100 fApplySPDDeadPbPb2011(kFALSE),
101 fApplySPDMisalignedPP2012(kFALSE),
102 fMaxDiffTRKV0Centr(-1.),
103 fRemoveTrackletOutliers(kFALSE),
106 fUseTrackSelectionWithFilterBits(kTRUE),
107 fUseCentrFlatteningInMC(kFALSE),
108 fHistCentrDistr(0x0),
109 fCutRatioClsOverCrossRowsTPC(0),
110 fCutRatioSignalNOverCrossRowsTPC(0),
111 fCutMinCrossedRowsTPCPtDep(""),
112 f1CutMinNCrossedRowsTPCPtDep(0x0)
115 // Default Constructor
117 fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
119 //--------------------------------------------------------------------------
120 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
121 AliAnalysisCuts(source),
122 fMinVtxType(source.fMinVtxType),
123 fMinVtxContr(source.fMinVtxContr),
124 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
125 fMaxVtxZ(source.fMaxVtxZ),
126 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
127 fTriggerMask(source.fTriggerMask),
128 fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
131 fnPtBins(source.fnPtBins),
132 fnPtBinLimits(source.fnPtBinLimits),
134 fnVars(source.fnVars),
136 fnVarsForOpt(source.fnVarsForOpt),
138 fGlobalIndex(source.fGlobalIndex),
141 fUsePID(source.fUsePID),
142 fUseAOD049(source.fUseAOD049),
144 fWhyRejection(source.fWhyRejection),
145 fEvRejectionBits(source.fEvRejectionBits),
146 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
147 fRecomputePrimVertex(source.fRecomputePrimVertex),
148 fUseMCVertex(source.fUseMCVertex),
149 fUsePhysicsSelection(source.fUsePhysicsSelection),
150 fOptPileup(source.fOptPileup),
151 fMinContrPileup(source.fMinContrPileup),
152 fMinDzPileup(source.fMinDzPileup),
153 fUseCentrality(source.fUseCentrality),
154 fMinCentrality(source.fMinCentrality),
155 fMaxCentrality(source.fMaxCentrality),
156 fFixRefs(source.fFixRefs),
157 fIsSelectedCuts(source.fIsSelectedCuts),
158 fIsSelectedPID(source.fIsSelectedPID),
159 fMinPtCand(source.fMinPtCand),
160 fMaxPtCand(source.fMaxPtCand),
161 fMaxRapidityCand(source.fMaxRapidityCand),
162 fKeepSignalMC(source.fKeepSignalMC),
163 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
164 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
165 fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
166 fApplySPDMisalignedPP2012(source.fApplySPDMisalignedPP2012),
167 fMaxDiffTRKV0Centr(source.fMaxDiffTRKV0Centr),
168 fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
169 fCutOnzVertexSPD(source.fCutOnzVertexSPD),
170 fKinkReject(source.fKinkReject),
171 fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
172 fUseCentrFlatteningInMC(source.fUseCentrFlatteningInMC),
173 fHistCentrDistr(0x0),
174 fCutRatioClsOverCrossRowsTPC(source.fCutRatioClsOverCrossRowsTPC),
175 fCutRatioSignalNOverCrossRowsTPC(source.fCutRatioSignalNOverCrossRowsTPC),
176 fCutMinCrossedRowsTPCPtDep(""),
177 f1CutMinNCrossedRowsTPCPtDep(0x0)
182 cout<<"Copy constructor"<<endl;
183 fTriggerClass[0] = source.fTriggerClass[0];
184 fTriggerClass[1] = source.fTriggerClass[1];
185 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
186 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
187 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
188 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
189 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
190 if(source.fPidHF) SetPidHF(source.fPidHF);
191 if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
192 if(source.fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
193 if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
197 //--------------------------------------------------------------------------
198 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
201 // assignment operator
203 if(&source == this) return *this;
205 AliAnalysisCuts::operator=(source);
207 fMinVtxType=source.fMinVtxType;
208 fMinVtxContr=source.fMinVtxContr;
209 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
210 fMaxVtxZ=source.fMaxVtxZ;
211 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
212 fTriggerMask=source.fTriggerMask;
213 fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
214 fTriggerClass[0]=source.fTriggerClass[0];
215 fTriggerClass[1]=source.fTriggerClass[1];
216 fnPtBins=source.fnPtBins;
217 fnPtBinLimits=source.fnPtBinLimits;
218 fnVars=source.fnVars;
219 fGlobalIndex=source.fGlobalIndex;
220 fnVarsForOpt=source.fnVarsForOpt;
221 fUsePID=source.fUsePID;
222 fUseAOD049=source.fUseAOD049;
223 if(fPidHF) delete fPidHF;
224 fPidHF=new AliAODPidHF(*(source.GetPidHF()));
225 fWhyRejection=source.fWhyRejection;
226 fEvRejectionBits=source.fEvRejectionBits;
227 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
228 fRecomputePrimVertex=source.fRecomputePrimVertex;
229 fUseMCVertex=source.fUseMCVertex;
230 fUsePhysicsSelection=source.fUsePhysicsSelection;
231 fOptPileup=source.fOptPileup;
232 fMinContrPileup=source.fMinContrPileup;
233 fMinDzPileup=source.fMinDzPileup;
234 fUseCentrality=source.fUseCentrality;
235 fMinCentrality=source.fMinCentrality;
236 fMaxCentrality=source.fMaxCentrality;
237 fFixRefs=source.fFixRefs;
238 fIsSelectedCuts=source.fIsSelectedCuts;
239 fIsSelectedPID=source.fIsSelectedPID;
240 fMinPtCand=source.fMinPtCand;
241 fMaxPtCand=source.fMaxPtCand;
242 fMaxRapidityCand=source.fMaxRapidityCand;
243 fKeepSignalMC=source.fKeepSignalMC;
244 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
245 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
246 fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
247 fApplySPDMisalignedPP2012=source.fApplySPDMisalignedPP2012;
248 fMaxDiffTRKV0Centr=source.fMaxDiffTRKV0Centr;
249 fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
250 fCutOnzVertexSPD=source.fCutOnzVertexSPD;
251 fKinkReject=source.fKinkReject;
252 fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
253 if(fHistCentrDistr) delete fHistCentrDistr;
254 fUseCentrFlatteningInMC=source.fUseCentrFlatteningInMC;
255 if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
257 if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
258 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
259 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
260 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
261 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
263 if(fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
264 if(f1CutMinNCrossedRowsTPCPtDep) delete f1CutMinNCrossedRowsTPCPtDep;
265 if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
266 fCutRatioClsOverCrossRowsTPC=source.fCutRatioClsOverCrossRowsTPC;
267 fCutRatioSignalNOverCrossRowsTPC=source.fCutRatioSignalNOverCrossRowsTPC;
272 //--------------------------------------------------------------------------
273 AliRDHFCuts::~AliRDHFCuts() {
275 // Default Destructor
277 if(fTrackCuts) { delete fTrackCuts; fTrackCuts=0; }
278 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
279 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
280 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
285 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
290 if(fHistCentrDistr)delete fHistCentrDistr;
292 if(f1CutMinNCrossedRowsTPCPtDep) {
293 delete f1CutMinNCrossedRowsTPCPtDep;
294 f1CutMinNCrossedRowsTPCPtDep = 0;
298 //---------------------------------------------------------------------------
299 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
301 // Centrality selection
303 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
304 AliWarning("Centrality estimator not valid");
307 Float_t centvalue=GetCentrality((AliAODEvent*)event);
308 if (centvalue<-998.){//-999 if no centralityP
311 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
314 // selection to flatten the centrality distribution
316 if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
324 //-------------------------------------------------
325 void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
326 // set the histo for centrality flattening
327 // the centrality is flatten in the range minCentr,maxCentr
328 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
329 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
330 // 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)
331 // 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
333 if(maxCentr<minCentr){
334 AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
337 if(fHistCentrDistr)delete fHistCentrDistr;
338 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
339 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
340 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
341 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
342 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
343 Double_t ref=0.,bincont=0.,binrefwidth=1.;
345 if(TMath::Abs(centrRef)<0.0001){
346 binref=fHistCentrDistr->GetMinimumBin();
347 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
348 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
350 else if(centrRef>0.){
351 binref=h->FindBin(centrRef);
352 if(binref<1||binref>h->GetNbinsX()){
353 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
355 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
356 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
359 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
360 binref=fHistCentrDistr->GetMaximumBin();
361 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
362 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
365 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
366 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
367 bincont=h->GetBinContent(j);
368 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
369 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
372 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
376 fHistCentrDistr->SetBinContent(0,switchTRand);
381 //-------------------------------------------------
382 Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
384 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
385 // Can be faster if it was required that fHistCentrDistr covers
386 // exactly the desired centrality range (e.g. part of the lines below should be done during the
387 // setting of the histo) and TH1::SetMinimum called
390 if(!fHistCentrDistr) return kTRUE;
391 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
392 // if(maxbin>fHistCentrDistr->GetNbinsX()){
393 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
396 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
397 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
398 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
400 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
401 if(gRandom->Uniform(1.)<bincont)return kTRUE;
405 if(centDigits*100.<bincont)return kTRUE;
409 //---------------------------------------------------------------------------
410 void AliRDHFCuts::SetupPID(AliVEvent *event) {
411 // Set the PID response object in the AliAODPidHF
412 // in case of old PID sets the TPC dE/dx BB parameterization
415 if(fPidHF->GetPidResponse()==0x0){
416 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
417 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
418 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
419 fPidHF->SetPidResponse(pidResp);
421 if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
422 if(fPidHF->GetOldPid()) {
425 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
426 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
428 // pp, from LHC10d onwards
429 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
430 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
431 // pp, 2011 low energy run
432 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
433 fPidHF->SetppLowEn2011(kTRUE);
434 fPidHF->SetOnePad(kFALSE);
437 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
439 if(isMC) fPidHF->SetMC(kTRUE);
440 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
441 fPidHF->SetMClowenpp2011(kTRUE);
442 fPidHF->SetBetheBloch();
444 // check that AliPIDResponse object was properly set in case of using OADB
445 if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
449 //---------------------------------------------------------------------------
450 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
454 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
462 if(fRecomputePrimVertex){
463 Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
472 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
473 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
479 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
480 // don't do for MC and for PbPb 2010 data
481 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
482 if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
483 (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
485 fEvRejectionBits+=1<<kNotSelTrigger;
490 // TEMPORARY FIX FOR GetEvent
491 Int_t nTracks=((AliAODEvent*)event)->GetNumberOfTracks();
492 for(Int_t itr=0; itr<nTracks; itr++){
493 AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
494 tr->SetAODEvent((AliAODEvent*)event);
497 // TEMPORARY FIX FOR REFERENCES
498 // Fix references to daughter tracks
500 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
501 // fixer->FixReferences((AliAODEvent*)event);
507 // physics selection requirements
508 if(fUsePhysicsSelection){
509 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
511 if(accept) fWhyRejection=7;
512 fEvRejectionBits+=1<<kPhysicsSelection;
515 if(fUseOnlyOneTrigger){
516 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
517 if(accept) fWhyRejection=7;
518 fEvRejectionBits+=1<<kPhysicsSelection;
525 // vertex requirements
527 const AliVVertex *vertex = event->GetPrimaryVertex();
531 fEvRejectionBits+=1<<kNoVertex;
533 TString title=vertex->GetTitle();
534 if(title.Contains("Z") && fMinVtxType>1){
536 fEvRejectionBits+=1<<kNoVertex;
538 else if(title.Contains("3D") && fMinVtxType>2){
540 fEvRejectionBits+=1<<kNoVertex;
542 if(vertex->GetNContributors()<fMinVtxContr){
544 fEvRejectionBits+=1<<kTooFewVtxContrib;
546 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
547 fEvRejectionBits+=1<<kZVtxOutFid;
548 if(accept) fWhyRejection=6;
553 if(fCutOnzVertexSPD>0){
554 const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
555 if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
557 fEvRejectionBits+=1<<kBadSPDVertex;
559 if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
560 fEvRejectionBits+=1<<kZVtxSPDOutFid;
561 if(accept) fWhyRejection=6;
564 if(fCutOnzVertexSPD==2 && vertex){
565 if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
566 fEvRejectionBits+=1<<kZVtxSPDOutFid;
567 if(accept) fWhyRejection=6;
575 if(fOptPileup==kRejectPileupEvent){
576 Int_t cutc=(Int_t)fMinContrPileup;
577 Double_t cutz=(Double_t)fMinDzPileup;
578 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
579 if(accept) fWhyRejection=1;
580 fEvRejectionBits+=1<<kPileup;
584 else if(fOptPileup==kRejectMVPileupEvent){
585 AliAnalysisUtils utils;
586 Bool_t isPUMV = utils.IsPileUpMV(event);
588 if(accept) fWhyRejection=1;
589 fEvRejectionBits+=1<<kPileup;
594 // centrality selection
595 if (fUseCentrality!=kCentOff) {
596 Int_t rejection=IsEventSelectedInCentrality(event);
597 Bool_t okCent=kFALSE;
598 if(rejection==0) okCent=kTRUE;
599 if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
601 if(accept) fWhyRejection=rejection;
602 if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
603 else fEvRejectionBits+=1<<kOutsideCentrality;
609 // PbPb2011 outliers in tracklets vs. VZERO and centTRK vs. centV0
610 if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
611 if(fRemoveTrackletOutliers){
612 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
613 Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
614 Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
615 if(ntracklets<1000. && v0cent<cutval){
616 if(accept) fWhyRejection=2;
617 fEvRejectionBits+=1<<kOutsideCentrality;
621 if(fMaxDiffTRKV0Centr>0.){
622 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
623 Double_t trkcent=GetCentrality((AliAODEvent*)event,kCentTRK);
624 if(TMath::Abs(trkcent-v0cent)>fMaxDiffTRKV0Centr){
625 if(accept) fWhyRejection=1;
626 fEvRejectionBits+=1<<kBadTrackV0Correl;
632 // Correcting PP2012 flag to remoce tracks crossing SPD misaligned staves for periods 12def
633 if(fApplySPDMisalignedPP2012 && !(event->GetRunNumber()>=195681 && event->GetRunNumber()<=197388)) fApplySPDMisalignedPP2012=false;
637 //---------------------------------------------------------------------------
638 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const{
640 // Daughter tracks selection
642 if(!fTrackCuts) return kTRUE;
644 Int_t ndaughters = d->GetNDaughters();
645 AliAODVertex *vAOD = d->GetPrimaryVtx();
646 Double_t pos[3],cov[6];
648 vAOD->GetCovarianceMatrix(cov);
649 const AliESDVertex vESD(pos,cov,100.,100);
653 for(Int_t idg=0; idg<ndaughters; idg++) {
654 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
655 if(!dgTrack) {retval = kFALSE; continue;}
656 //printf("charge %d\n",dgTrack->Charge());
657 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
659 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
660 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
662 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
667 //---------------------------------------------------------------------------
668 Bool_t AliRDHFCuts::CheckPtDepCrossedRows(TString rows,Bool_t print) const {
670 // Check the correctness of the string syntax
674 if(!rows.Contains("pt")) {
675 if(print) AliError("string must contain \"pt\"");
680 //---------------------------------------------------------------------------
681 void AliRDHFCuts::SetMinCrossedRowsTPCPtDep(const char *rows){
683 //Create the TFormula from TString for TPC crossed rows pT dependent cut
687 // setting data member that describes the TPC crossed rows pT dependent cut
688 fCutMinCrossedRowsTPCPtDep = rows;
690 // creating TFormula from TString
691 if(f1CutMinNCrossedRowsTPCPtDep){
692 delete f1CutMinNCrossedRowsTPCPtDep;
693 // resetting TFormula
694 f1CutMinNCrossedRowsTPCPtDep = 0;
696 if(!CheckPtDepCrossedRows(rows,kTRUE))return;
699 tmp.ReplaceAll("pt","x");
700 f1CutMinNCrossedRowsTPCPtDep = new TFormula("f1CutMinNCrossedRowsTPCPtDep",tmp.Data());
704 //---------------------------------------------------------------------------
705 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const{
707 // Convert to ESDtrack, relate to vertex and check cuts
709 if(!cuts) return kTRUE;
711 if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
714 // convert to ESD track here
715 AliESDtrack esdTrack(track);
716 // set the TPC cluster info
717 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
718 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
719 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
720 // needed to calculate the impact parameters
721 esdTrack.RelateToVertex(primary,0.,3.);
723 //applying ESDtrackCut
724 if(!cuts->IsSelected(&esdTrack)) return kFALSE;
726 //appliyng kink rejection
728 AliAODVertex *maybeKink=track->GetProdVertex();
729 if(maybeKink->GetType()==AliAODVertex::kKink) return kFALSE;
732 //appliyng TPC crossed rows pT dependent cut
733 if(f1CutMinNCrossedRowsTPCPtDep){
734 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
735 if(nCrossedRowsTPC<f1CutMinNCrossedRowsTPCPtDep->Eval(esdTrack.Pt())) return kFALSE;
738 //appliyng NTPCcls/NTPCcrossedRows cut
739 if(fCutRatioClsOverCrossRowsTPC){
740 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
741 Float_t nClustersTPC = esdTrack.GetTPCNcls();
742 if(nCrossedRowsTPC!=0){
743 Float_t ratio = nClustersTPC/nCrossedRowsTPC;
744 if(ratio<fCutRatioClsOverCrossRowsTPC) return kFALSE;
749 //appliyng TPCsignalN/NTPCcrossedRows cut
750 if(fCutRatioSignalNOverCrossRowsTPC){
751 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
752 Float_t nTPCsignal = esdTrack.GetTPCsignalN();
753 if(nCrossedRowsTPC!=0){
754 Float_t ratio = nTPCsignal/nCrossedRowsTPC;
755 if(ratio<fCutRatioSignalNOverCrossRowsTPC) return kFALSE;
761 if(fOptPileup==kRejectTracksFromPileupVertex){
763 // we need either to have here the AOD Event,
764 // or to have the pileup vertex object
767 if(fApplySPDDeadPbPb2011){
768 Bool_t deadSPDLay1PbPb2011[20][4]={
769 {kTRUE,kTRUE,kTRUE,kTRUE},
770 {kTRUE,kTRUE,kTRUE,kTRUE},
771 {kTRUE,kTRUE,kTRUE,kTRUE},
772 {kTRUE,kTRUE,kTRUE,kTRUE},
773 {kTRUE,kTRUE,kTRUE,kTRUE},
774 {kFALSE,kFALSE,kTRUE,kTRUE},
775 {kTRUE,kTRUE,kFALSE,kFALSE},
776 {kTRUE,kTRUE,kTRUE,kTRUE},
777 {kFALSE,kFALSE,kTRUE,kTRUE},
778 {kTRUE,kTRUE,kTRUE,kTRUE},
779 {kTRUE,kTRUE,kFALSE,kFALSE},
780 {kTRUE,kTRUE,kTRUE,kTRUE},
781 {kFALSE,kFALSE,kFALSE,kFALSE},
782 {kFALSE,kFALSE,kTRUE,kTRUE},
783 {kFALSE,kFALSE,kFALSE,kFALSE},
784 {kFALSE,kFALSE,kFALSE,kFALSE},
785 {kTRUE,kTRUE,kTRUE,kTRUE},
786 {kTRUE,kTRUE,kFALSE,kFALSE},
787 {kFALSE,kFALSE,kFALSE,kFALSE},
788 {kFALSE,kFALSE,kFALSE,kFALSE}
790 Bool_t deadSPDLay2PbPb2011[40][4]={
791 {kTRUE,kTRUE,kTRUE,kTRUE},
792 {kTRUE,kTRUE,kTRUE,kTRUE},
793 {kTRUE,kTRUE,kTRUE,kTRUE},
794 {kTRUE,kTRUE,kTRUE,kTRUE},
795 {kTRUE,kTRUE,kTRUE,kTRUE},
796 {kTRUE,kTRUE,kTRUE,kTRUE},
797 {kTRUE,kTRUE,kTRUE,kTRUE},
798 {kTRUE,kTRUE,kTRUE,kTRUE},
799 {kTRUE,kTRUE,kTRUE,kTRUE},
800 {kTRUE,kTRUE,kTRUE,kTRUE},
801 {kTRUE,kTRUE,kTRUE,kTRUE},
802 {kTRUE,kTRUE,kTRUE,kTRUE},
803 {kFALSE,kFALSE,kFALSE,kFALSE},
804 {kFALSE,kFALSE,kTRUE,kTRUE},
805 {kTRUE,kTRUE,kTRUE,kTRUE},
806 {kTRUE,kTRUE,kTRUE,kTRUE},
807 {kTRUE,kTRUE,kFALSE,kFALSE},
808 {kTRUE,kTRUE,kTRUE,kTRUE},
809 {kTRUE,kTRUE,kTRUE,kTRUE},
810 {kTRUE,kTRUE,kTRUE,kTRUE},
811 {kFALSE,kFALSE,kFALSE,kFALSE},
812 {kFALSE,kFALSE,kFALSE,kFALSE},
813 {kTRUE,kTRUE,kTRUE,kTRUE},
814 {kTRUE,kTRUE,kTRUE,kTRUE},
815 {kFALSE,kFALSE,kFALSE,kFALSE},
816 {kFALSE,kFALSE,kFALSE,kFALSE},
817 {kTRUE,kTRUE,kTRUE,kTRUE},
818 {kTRUE,kTRUE,kTRUE,kTRUE},
819 {kFALSE,kFALSE,kFALSE,kFALSE},
820 {kFALSE,kFALSE,kFALSE,kFALSE},
821 {kFALSE,kFALSE,kFALSE,kFALSE},
822 {kFALSE,kFALSE,kFALSE,kFALSE},
823 {kTRUE,kTRUE,kTRUE,kTRUE},
824 {kTRUE,kTRUE,kTRUE,kTRUE},
825 {kTRUE,kTRUE,kFALSE,kFALSE},
826 {kTRUE,kTRUE,kTRUE,kTRUE},
827 {kFALSE,kFALSE,kFALSE,kFALSE},
828 {kFALSE,kFALSE,kFALSE,kFALSE},
829 {kFALSE,kFALSE,kFALSE,kFALSE},
830 {kFALSE,kFALSE,kFALSE,kFALSE}
832 Double_t xyz1[3],xyz2[3];
833 esdTrack.GetXYZAt(3.9,0.,xyz1);
834 esdTrack.GetXYZAt(7.6,0.,xyz2);
835 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
836 if(phi1<0) phi1+=2*TMath::Pi();
837 Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
838 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
839 if(phi2<0) phi2+=2*TMath::Pi();
840 Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
841 Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
842 Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
843 Bool_t lay1ok=kFALSE;
844 if(mod1>=0 && mod1<4 && lad1<20){
845 lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
847 Bool_t lay2ok=kFALSE;
848 if(mod2>=0 && mod2<4 && lad2<40){
849 lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
851 if(!lay1ok && !lay2ok) return kFALSE;
854 if(fApplySPDMisalignedPP2012) {
855 // Cut tracks crossing the SPD at 5.6<phi<2pi
856 Double_t xyz1[3],xyz2[3];
857 esdTrack.GetXYZAt(3.9,0.,xyz1);
858 esdTrack.GetXYZAt(7.6,0.,xyz2);
859 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
860 if(phi1<0) phi1+=2*TMath::Pi();
861 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
862 if(phi2<0) phi2+=2*TMath::Pi();
864 if(phi1>5.6 && phi1<2.*TMath::Pi()) lay1ok=kFALSE;
866 if(phi2>5.6 && phi2<2.*TMath::Pi()) lay2ok=kFALSE;
867 if(!lay1ok || !lay2ok) return kFALSE;
872 //---------------------------------------------------------------------------
873 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
877 delete [] fPtBinLimits;
879 printf("Changing the pt bins\n");
882 if(nPtBinLimits != fnPtBins+1){
883 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
884 SetNPtBins(nPtBinLimits-1);
887 fnPtBinLimits = nPtBinLimits;
889 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
890 fPtBinLimits = new Float_t[fnPtBinLimits];
891 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
895 //---------------------------------------------------------------------------
896 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
897 // Set the variable names
902 //printf("Changing the variable names\n");
905 printf("Wrong number of variables: it has to be %d\n",fnVars);
909 fVarNames = new TString[nVars];
910 fIsUpperCut = new Bool_t[nVars];
911 for(Int_t iv=0; iv<nVars; iv++) {
912 fVarNames[iv] = varNames[iv];
913 fIsUpperCut[iv] = isUpperCut[iv];
918 //---------------------------------------------------------------------------
919 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
920 // Set the variables to be used for cuts optimization
923 delete [] fVarsForOpt;
925 //printf("Changing the variables for cut optimization\n");
928 if(nVars==0){//!=fnVars) {
929 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
934 fVarsForOpt = new Bool_t[fnVars];
935 for(Int_t iv=0; iv<fnVars; iv++) {
936 fVarsForOpt[iv]=forOpt[iv];
937 if(fVarsForOpt[iv]) fnVarsForOpt++;
943 //---------------------------------------------------------------------------
944 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
946 // set centrality estimator
949 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
955 //---------------------------------------------------------------------------
956 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
961 printf("Wrong number of variables: it has to be %d\n",fnVars);
964 if(nPtBins!=fnPtBins) {
965 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
969 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
972 for(Int_t iv=0; iv<fnVars; iv++) {
974 for(Int_t ib=0; ib<fnPtBins; ib++) {
977 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
978 cout<<"Overflow, exit..."<<endl;
982 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
988 //---------------------------------------------------------------------------
989 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
993 if(glIndex != fGlobalIndex){
994 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
997 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
999 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
1000 fCutsRD[iGl] = cutsRDGlob[iGl];
1004 //---------------------------------------------------------------------------
1005 void AliRDHFCuts::PrintAll() const {
1007 // print all cuts values
1010 printf("Minimum vtx type %d\n",fMinVtxType);
1011 printf("Minimum vtx contr %d\n",fMinVtxContr);
1012 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
1013 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
1014 printf("Use PID %d OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
1015 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
1016 printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
1017 printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
1018 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
1019 if(fOptPileup==1) printf(" -- Reject pileup event");
1020 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
1021 if(fUseCentrality>0) {
1022 TString estimator="";
1023 if(fUseCentrality==1) estimator = "V0";
1024 if(fUseCentrality==2) estimator = "Tracks";
1025 if(fUseCentrality==3) estimator = "Tracklets";
1026 if(fUseCentrality==4) estimator = "SPD clusters outer";
1027 if(fUseCentrality==5) estimator = "ZNA";
1028 if(fUseCentrality==6) estimator = "ZPA";
1029 if(fUseCentrality==7) estimator = "V0A";
1030 printf("Centrality class considered: %.1f-%.1f, estimated with %s\n",fMinCentrality,fMaxCentrality,estimator.Data());
1032 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
1034 if(fCutRatioClsOverCrossRowsTPC) printf("N TPC Clusters > %f N TPC Crossed Rows\n", fCutRatioClsOverCrossRowsTPC);
1035 if(fCutRatioSignalNOverCrossRowsTPC) printf("N TPC Points for dE/dx > %f N TPC Crossed Rows\n", fCutRatioSignalNOverCrossRowsTPC);
1036 if(f1CutMinNCrossedRowsTPCPtDep) printf("N TPC Crossed Rows pT-dependent cut: %s\n", fCutMinCrossedRowsTPCPtDep.Data());
1039 cout<<"Array of variables"<<endl;
1040 for(Int_t iv=0;iv<fnVars;iv++){
1041 cout<<fVarNames[iv]<<"\t";
1046 cout<<"Array of optimization"<<endl;
1047 for(Int_t iv=0;iv<fnVars;iv++){
1048 cout<<fVarsForOpt[iv]<<"\t";
1053 cout<<"Array of upper/lower cut"<<endl;
1054 for(Int_t iv=0;iv<fnVars;iv++){
1055 cout<<fIsUpperCut[iv]<<"\t";
1060 cout<<"Array of ptbin limits"<<endl;
1061 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
1062 cout<<fPtBinLimits[ib]<<"\t";
1067 cout<<"Matrix of cuts"<<endl;
1068 for(Int_t iv=0;iv<fnVars;iv++){
1069 for(Int_t ib=0;ib<fnPtBins;ib++){
1070 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
1076 if(fPidHF) fPidHF->PrintAll();
1080 //--------------------------------------------------------------------------
1081 void AliRDHFCuts::PrintTrigger() const{
1082 // print the trigger selection
1084 printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
1086 cout<<" Trigger selection pattern: ";
1087 if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
1088 if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
1089 if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
1090 if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
1091 if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
1092 if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
1093 if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
1094 if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
1095 if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
1096 if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
1097 cout << endl<< endl;
1101 //---------------------------------------------------------------------------
1102 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
1107 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
1112 //cout<<"Initialization..."<<endl;
1113 cutsRD=new Float_t*[fnVars];
1114 for(iv=0; iv<fnVars; iv++) {
1115 cutsRD[iv] = new Float_t[fnPtBins];
1119 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
1120 GetVarPtIndex(iGlobal,iv,ib);
1121 cutsRD[iv][ib] = fCutsRD[iGlobal];
1127 //---------------------------------------------------------------------------
1128 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
1130 // give the global index from variable and pt bin
1132 return iPtBin*fnVars+iVar;
1135 //---------------------------------------------------------------------------
1136 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
1138 //give the index of the variable and of the pt bin from the global index
1140 iPtBin=(Int_t)iGlob/fnVars;
1146 //---------------------------------------------------------------------------
1147 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1149 //give the pt bin where the pt lies.
1152 if(pt<fPtBinLimits[0])return ptbin;
1153 for (Int_t i=0;i<fnPtBins;i++){
1154 if(pt<fPtBinLimits[i+1]) {
1161 //-------------------------------------------------------------------
1162 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1164 // Give the value of cut set for the variable iVar and the pt bin iPtBin
1167 cout<<"Cuts not iniziaisez yet"<<endl;
1170 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1172 //-------------------------------------------------------------------
1173 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1175 // Get centrality percentile
1178 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1179 if(mcArray) {fUseAOD049=kFALSE;}
1181 AliAODHeader *header=aodEvent->GetHeader();
1182 AliCentrality *centrality=header->GetCentralityP();
1184 Bool_t isSelRun=kFALSE;
1185 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1186 if(!centrality) return cent;
1188 if (estimator==kCentV0M){
1189 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1191 Int_t quality = centrality->GetQuality();
1192 if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1193 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1195 Int_t runnum=aodEvent->GetRunNumber();
1196 for(Int_t ir=0;ir<5;ir++){
1197 if(runnum==selRun[ir]){
1202 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1206 //temporary fix for AOD049 outliers
1207 if(fUseAOD049&¢>=0){
1209 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1210 v0+=aodV0->GetMTotV0A();
1211 v0+=aodV0->GetMTotV0C();
1212 if(cent==0&&v0<19500)return -1;//filtering issue
1213 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1214 Float_t val= 1.30552 + 0.147931 * v0;
1215 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};
1216 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1220 if (estimator==kCentTRK) {
1221 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1223 Int_t quality = centrality->GetQuality();
1225 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1227 Int_t runnum=aodEvent->GetRunNumber();
1228 for(Int_t ir=0;ir<5;ir++){
1229 if(runnum==selRun[ir]){
1234 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1239 if (estimator==kCentTKL){
1240 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1242 Int_t quality = centrality->GetQuality();
1244 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1246 Int_t runnum=aodEvent->GetRunNumber();
1247 for(Int_t ir=0;ir<5;ir++){
1248 if(runnum==selRun[ir]){
1253 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1258 if (estimator==kCentCL1){
1259 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1261 Int_t quality = centrality->GetQuality();
1263 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1265 Int_t runnum=aodEvent->GetRunNumber();
1266 for(Int_t ir=0;ir<5;ir++){
1267 if(runnum==selRun[ir]){
1272 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1277 if (estimator==kCentZNA){
1278 cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1280 Int_t quality = centrality->GetQuality();
1282 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1284 Int_t runnum=aodEvent->GetRunNumber();
1285 for(Int_t ir=0;ir<5;ir++){
1286 if(runnum==selRun[ir]){
1291 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1296 if (estimator==kCentZPA){
1297 cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1299 Int_t quality = centrality->GetQuality();
1301 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1303 Int_t runnum=aodEvent->GetRunNumber();
1304 for(Int_t ir=0;ir<5;ir++){
1305 if(runnum==selRun[ir]){
1310 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1315 if (estimator==kCentV0A){
1316 cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1318 Int_t quality = centrality->GetQuality();
1320 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1322 Int_t runnum=aodEvent->GetRunNumber();
1323 for(Int_t ir=0;ir<5;ir++){
1324 if(runnum==selRun[ir]){
1329 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1334 AliWarning("Centrality estimator not valid");
1346 //-------------------------------------------------------------------
1347 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1349 // Compare two cuts objects
1352 Bool_t areEqual=kTRUE;
1354 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1356 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1358 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1360 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1362 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1364 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1366 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1368 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1370 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1372 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;}
1376 for(Int_t iv=0;iv<fnVars;iv++) {
1377 for(Int_t ib=0;ib<fnPtBins;ib++) {
1378 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1379 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1388 //---------------------------------------------------------------------------
1389 void AliRDHFCuts::MakeTable() const {
1391 // print cuts values in table format
1394 TString ptString = "pT range";
1395 if(fVarNames && fPtBinLimits && fCutsRD){
1396 TString firstLine(Form("* %-15s",ptString.Data()));
1397 for (Int_t ivar=0; ivar<fnVars; ivar++){
1398 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1399 if (ivar == fnVars){
1403 Printf("%s",firstLine.Data());
1405 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1407 if (ipt==fnPtBins-1){
1408 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1411 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1413 for (Int_t ivar=0; ivar<fnVars; ivar++){
1414 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1416 Printf("%s",line.Data());
1424 //--------------------------------------------------------------------------
1425 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1426 AliAODEvent *aod) const
1429 // Recalculate primary vertex without daughters
1433 AliError("Can not remove daughters from vertex without AOD event");
1437 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1439 AliDebug(2,"Removal of daughter tracks failed");
1444 //set recalculed primary vertex
1445 d->SetOwnPrimaryVtx(recvtx);
1450 //--------------------------------------------------------------------------
1451 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1454 // Recalculate primary vertex without daughters
1458 AliError("Can not get MC vertex without AOD event");
1463 AliAODMCHeader *mcHeader =
1464 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1466 AliError("Can not get MC vertex without AODMCHeader event");
1470 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1471 mcHeader->GetVertex(pos);
1472 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1475 AliDebug(2,"Removal of daughter tracks failed");
1479 //set recalculed primary vertex
1480 d->SetOwnPrimaryVtx(recvtx);
1482 d->RecalculateImpPars(recvtx,aod);
1488 //--------------------------------------------------------------------------
1489 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1491 AliAODVertex *origownvtx) const
1494 // Clean-up own primary vertex if needed
1497 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1498 d->UnsetOwnPrimaryVtx();
1500 d->SetOwnPrimaryVtx(origownvtx);
1501 delete origownvtx; origownvtx=NULL;
1503 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1506 delete origownvtx; origownvtx=NULL;
1511 //--------------------------------------------------------------------------
1512 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1515 // Checks if this candidate is matched to MC signal
1518 if(!aod) return kFALSE;
1521 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1523 if(!mcArray) return kFALSE;
1526 Int_t label = d->MatchToMC(pdg,mcArray);
1529 //printf("MATCH!\n");
1537 //--------------------------------------------------------------------------
1538 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1539 // recompute event primary vertex from AOD tracks
1541 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1542 vertexer->SetITSMode();
1543 vertexer->SetMinClusters(3);
1545 AliAODVertex* pvtx=event->GetPrimaryVertex();
1546 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1547 Float_t diamondcovxy[3];
1548 event->GetDiamondCovXY(diamondcovxy);
1549 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1550 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1551 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1552 vertexer->SetVtxStart(diamond);
1553 delete diamond; diamond=NULL;
1556 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1557 if(!vertexESD) return kFALSE;
1558 if(vertexESD->GetNContributors()<=0) {
1559 //AliDebug(2,"vertexing failed");
1560 delete vertexESD; vertexESD=NULL;
1563 delete vertexer; vertexer=NULL;
1565 // convert to AliAODVertex
1566 Double_t pos[3],cov[6],chi2perNDF;
1567 vertexESD->GetXYZ(pos); // position
1568 vertexESD->GetCovMatrix(cov); //covariance matrix
1569 chi2perNDF = vertexESD->GetChi2toNDF();
1570 delete vertexESD; vertexESD=NULL;
1572 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1573 pvtx->SetChi2perNDF(chi2perNDF);
1574 pvtx->SetCovMatrix(cov);