]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliRDHFCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCuts.cxx
CommitLineData
0c22e2ac 1/**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18/////////////////////////////////////////////////////////////
19//
20// Base class for cuts on AOD reconstructed heavy-flavour decay
21//
22// Author: A.Dainese, andrea.dainese@pd.infn.it
23/////////////////////////////////////////////////////////////
24#include <Riostream.h>
25
26#include "AliVEvent.h"
27#include "AliESDEvent.h"
28#include "AliAODEvent.h"
29#include "AliVVertex.h"
30#include "AliESDVertex.h"
31#include "AliLog.h"
32#include "AliAODVertex.h"
33#include "AliESDtrack.h"
34#include "AliAODTrack.h"
35#include "AliESDtrackCuts.h"
36#include "AliCentrality.h"
37#include "AliAODRecoDecayHF.h"
38#include "AliAnalysisVertexingHF.h"
39#include "AliAODMCHeader.h"
40#include "AliAODMCParticle.h"
4bdc0315 41#include "AliVertexerTracks.h"
0c22e2ac 42#include "AliRDHFCuts.h"
43#include "AliAnalysisManager.h"
44#include "AliInputEventHandler.h"
45#include "AliPIDResponse.h"
5d3cf93b 46#include "AliAnalysisUtils.h"
db0ff379 47#include "TRandom.h"
09bfd399 48#include <TF1.h>
0c22e2ac 49
c64cb1f6 50using std::cout;
51using std::endl;
52
0c22e2ac 53ClassImp(AliRDHFCuts)
54
55//--------------------------------------------------------------------------
56AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
57AliAnalysisCuts(name,title),
58fMinVtxType(3),
59fMinVtxContr(1),
60fMaxVtxRedChi2(1e6),
61fMaxVtxZ(10.),
62fMinSPDMultiplicity(0),
63fTriggerMask(AliVEvent::kAnyINT),
64fUseOnlyOneTrigger(kFALSE),
0c22e2ac 65fTrackCuts(0),
66fnPtBins(1),
67fnPtBinLimits(1),
68fPtBinLimits(0),
69fnVars(1),
70fVarNames(0),
71fnVarsForOpt(0),
72fVarsForOpt(0),
73fGlobalIndex(1),
74fCutsRD(0),
75fIsUpperCut(0),
76fUsePID(kFALSE),
77fUseAOD049(kFALSE),
78fPidHF(0),
79fWhyRejection(0),
80fEvRejectionBits(0),
81fRemoveDaughtersFromPrimary(kFALSE),
82fUseMCVertex(kFALSE),
83fUsePhysicsSelection(kTRUE),
84fOptPileup(0),
85fMinContrPileup(3),
86fMinDzPileup(0.6),
87fUseCentrality(0),
88fMinCentrality(0.),
89fMaxCentrality(100.),
90fFixRefs(kFALSE),
91fIsSelectedCuts(0),
92fIsSelectedPID(0),
93fMinPtCand(-1.),
94fMaxPtCand(100000.),
af5d687e 95fMaxRapidityCand(-999.),
0c22e2ac 96fKeepSignalMC(kFALSE),
97fIsCandTrackSPDFirst(kFALSE),
851d1cae 98fMaxPtCandTrackSPDFirst(0.),
99fApplySPDDeadPbPb2011(kFALSE),
180a04cb 100fApplySPDMisalignedPP2012(kFALSE),
8da55fd2 101fMaxDiffTRKV0Centr(-1.),
fdd5a897 102fRemoveTrackletOutliers(kFALSE),
4346ffd2 103fCutOnzVertexSPD(0),
3e075b37 104fKinkReject(kFALSE),
09bfd399 105fUseTrackSelectionWithFilterBits(kTRUE),
106fUseCentrFlatteningInMC(kFALSE),
107fHistCentrDistr(0x0),
108fCutRatioClsOverCrossRowsTPC(0),
109fCutRatioSignalNOverCrossRowsTPC(0),
110fCutMinCrossedRowsTPCPtDep(""),
111f1CutMinNCrossedRowsTPCPtDep(0x0)
0c22e2ac 112{
113 //
114 // Default Constructor
115 //
68c15f66 116 fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
0c22e2ac 117}
118//--------------------------------------------------------------------------
119AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
120 AliAnalysisCuts(source),
121 fMinVtxType(source.fMinVtxType),
122 fMinVtxContr(source.fMinVtxContr),
123 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
124 fMaxVtxZ(source.fMaxVtxZ),
125 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
126 fTriggerMask(source.fTriggerMask),
127 fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
68c15f66 128 fTriggerClass(),
0c22e2ac 129 fTrackCuts(0),
130 fnPtBins(source.fnPtBins),
131 fnPtBinLimits(source.fnPtBinLimits),
132 fPtBinLimits(0),
133 fnVars(source.fnVars),
134 fVarNames(0),
135 fnVarsForOpt(source.fnVarsForOpt),
136 fVarsForOpt(0),
137 fGlobalIndex(source.fGlobalIndex),
138 fCutsRD(0),
139 fIsUpperCut(0),
140 fUsePID(source.fUsePID),
141 fUseAOD049(source.fUseAOD049),
142 fPidHF(0),
143 fWhyRejection(source.fWhyRejection),
144 fEvRejectionBits(source.fEvRejectionBits),
145 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
146 fUseMCVertex(source.fUseMCVertex),
147 fUsePhysicsSelection(source.fUsePhysicsSelection),
148 fOptPileup(source.fOptPileup),
149 fMinContrPileup(source.fMinContrPileup),
150 fMinDzPileup(source.fMinDzPileup),
151 fUseCentrality(source.fUseCentrality),
152 fMinCentrality(source.fMinCentrality),
153 fMaxCentrality(source.fMaxCentrality),
154 fFixRefs(source.fFixRefs),
155 fIsSelectedCuts(source.fIsSelectedCuts),
156 fIsSelectedPID(source.fIsSelectedPID),
157 fMinPtCand(source.fMinPtCand),
158 fMaxPtCand(source.fMaxPtCand),
af5d687e 159 fMaxRapidityCand(source.fMaxRapidityCand),
0c22e2ac 160 fKeepSignalMC(source.fKeepSignalMC),
161 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
851d1cae 162 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
163 fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
180a04cb 164 fApplySPDMisalignedPP2012(source.fApplySPDMisalignedPP2012),
8da55fd2 165 fMaxDiffTRKV0Centr(source.fMaxDiffTRKV0Centr),
fdd5a897 166 fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
4346ffd2 167 fCutOnzVertexSPD(source.fCutOnzVertexSPD),
3e075b37 168 fKinkReject(source.fKinkReject),
db0ff379 169 fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
30399c33 170 fUseCentrFlatteningInMC(source.fUseCentrFlatteningInMC),
09bfd399 171 fHistCentrDistr(0x0),
172 fCutRatioClsOverCrossRowsTPC(source.fCutRatioClsOverCrossRowsTPC),
173 fCutRatioSignalNOverCrossRowsTPC(source.fCutRatioSignalNOverCrossRowsTPC),
849b4fa2 174 fCutMinCrossedRowsTPCPtDep(""),
09bfd399 175 f1CutMinNCrossedRowsTPCPtDep(0x0)
0c22e2ac 176{
177 //
178 // Copy constructor
179 //
180 cout<<"Copy constructor"<<endl;
68c15f66 181 fTriggerClass[0] = source.fTriggerClass[0];
182 fTriggerClass[1] = source.fTriggerClass[1];
0c22e2ac 183 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
184 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
185 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
186 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
187 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
188 if(source.fPidHF) SetPidHF(source.fPidHF);
db0ff379 189 if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
09bfd399 190 if(source.fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
191 if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
0c22e2ac 192 PrintAll();
193
194}
195//--------------------------------------------------------------------------
196AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
197{
198 //
199 // assignment operator
200 //
201 if(&source == this) return *this;
202
203 AliAnalysisCuts::operator=(source);
204
205 fMinVtxType=source.fMinVtxType;
206 fMinVtxContr=source.fMinVtxContr;
207 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
208 fMaxVtxZ=source.fMaxVtxZ;
209 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
210 fTriggerMask=source.fTriggerMask;
211 fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
68c15f66 212 fTriggerClass[0]=source.fTriggerClass[0];
213 fTriggerClass[1]=source.fTriggerClass[1];
0c22e2ac 214 fnPtBins=source.fnPtBins;
215 fnPtBinLimits=source.fnPtBinLimits;
216 fnVars=source.fnVars;
217 fGlobalIndex=source.fGlobalIndex;
218 fnVarsForOpt=source.fnVarsForOpt;
219 fUsePID=source.fUsePID;
220 fUseAOD049=source.fUseAOD049;
221 if(fPidHF) delete fPidHF;
222 fPidHF=new AliAODPidHF(*(source.GetPidHF()));
223 fWhyRejection=source.fWhyRejection;
224 fEvRejectionBits=source.fEvRejectionBits;
225 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
226 fUseMCVertex=source.fUseMCVertex;
227 fUsePhysicsSelection=source.fUsePhysicsSelection;
228 fOptPileup=source.fOptPileup;
229 fMinContrPileup=source.fMinContrPileup;
230 fMinDzPileup=source.fMinDzPileup;
231 fUseCentrality=source.fUseCentrality;
232 fMinCentrality=source.fMinCentrality;
233 fMaxCentrality=source.fMaxCentrality;
234 fFixRefs=source.fFixRefs;
235 fIsSelectedCuts=source.fIsSelectedCuts;
236 fIsSelectedPID=source.fIsSelectedPID;
237 fMinPtCand=source.fMinPtCand;
238 fMaxPtCand=source.fMaxPtCand;
af5d687e 239 fMaxRapidityCand=source.fMaxRapidityCand;
0c22e2ac 240 fKeepSignalMC=source.fKeepSignalMC;
241 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
242 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
851d1cae 243 fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
180a04cb 244 fApplySPDMisalignedPP2012=source.fApplySPDMisalignedPP2012;
8da55fd2 245 fMaxDiffTRKV0Centr=source.fMaxDiffTRKV0Centr;
851d1cae 246 fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
4346ffd2 247 fCutOnzVertexSPD=source.fCutOnzVertexSPD;
fdd5a897 248 fKinkReject=source.fKinkReject;
3e075b37 249 fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
db0ff379 250 if(fHistCentrDistr) delete fHistCentrDistr;
30399c33 251 fUseCentrFlatteningInMC=source.fUseCentrFlatteningInMC;
db0ff379 252 if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
0c22e2ac 253
254 if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
255 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
256 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
257 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
258 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
09bfd399 259
260 if(fCutMinCrossedRowsTPCPtDep) fCutMinCrossedRowsTPCPtDep=source.fCutMinCrossedRowsTPCPtDep;
261 if(f1CutMinNCrossedRowsTPCPtDep) delete f1CutMinNCrossedRowsTPCPtDep;
262 if(source.f1CutMinNCrossedRowsTPCPtDep) f1CutMinNCrossedRowsTPCPtDep=new TFormula(*(source.f1CutMinNCrossedRowsTPCPtDep));
263 fCutRatioClsOverCrossRowsTPC=source.fCutRatioClsOverCrossRowsTPC;
264 fCutRatioSignalNOverCrossRowsTPC=source.fCutRatioSignalNOverCrossRowsTPC;
0c22e2ac 265 PrintAll();
266
267 return *this;
268}
269//--------------------------------------------------------------------------
270AliRDHFCuts::~AliRDHFCuts() {
271 //
272 // Default Destructor
273 //
85204d46 274 if(fTrackCuts) { delete fTrackCuts; fTrackCuts=0; }
0c22e2ac 275 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
276 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
277 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
278 if(fCutsRD) {
279 delete [] fCutsRD;
280 fCutsRD=0;
281 }
282 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
283 if(fPidHF){
284 delete fPidHF;
285 fPidHF=0;
286 }
db0ff379 287 if(fHistCentrDistr)delete fHistCentrDistr;
09bfd399 288
289 if(f1CutMinNCrossedRowsTPCPtDep) {
290 delete f1CutMinNCrossedRowsTPCPtDep;
291 f1CutMinNCrossedRowsTPCPtDep = 0;
292 }
293
0c22e2ac 294}
295//---------------------------------------------------------------------------
296Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
297 //
298 // Centrality selection
299 //
300 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
301 AliWarning("Centrality estimator not valid");
302 return 3;
303 }else{
304 Float_t centvalue=GetCentrality((AliAODEvent*)event);
305 if (centvalue<-998.){//-999 if no centralityP
306 return 0;
307 }else{
308 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
309 return 2;
db0ff379 310 }
311 // selection to flatten the centrality distribution
312 if(fHistCentrDistr){
313 if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
314 }
0c22e2ac 315 }
316 }
317 return 0;
318}
db0ff379 319
320
321//-------------------------------------------------
322void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
323 // set the histo for centrality flattening
324 // the centrality is flatten in the range minCentr,maxCentr
325 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
326 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
327 // 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)
328 // 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
329
330 if(maxCentr<minCentr){
331 AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
332 }
333
334 if(fHistCentrDistr)delete fHistCentrDistr;
335 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
336 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
337 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
338 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
339 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
340 Double_t ref=0.,bincont=0.,binrefwidth=1.;
341 Int_t binref=0;
342 if(TMath::Abs(centrRef)<0.0001){
343 binref=fHistCentrDistr->GetMinimumBin();
344 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
345 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
346 }
347 else if(centrRef>0.){
348 binref=h->FindBin(centrRef);
349 if(binref<1||binref>h->GetNbinsX()){
350 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
351 }
352 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
353 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
354 }
355 else{
356 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
357 binref=fHistCentrDistr->GetMaximumBin();
358 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
359 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
360 }
361
362 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
363 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
364 bincont=h->GetBinContent(j);
365 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
366 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
367 }
368 else{
369 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
370 }
371 }
372
373 fHistCentrDistr->SetBinContent(0,switchTRand);
374 return;
375
376}
377
378//-------------------------------------------------
379Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
380 //
381 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
382 // Can be faster if it was required that fHistCentrDistr covers
383 // exactly the desired centrality range (e.g. part of the lines below should be done during the
384 // setting of the histo) and TH1::SetMinimum called
385 //
386
387 if(!fHistCentrDistr) return kTRUE;
388 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
389 // if(maxbin>fHistCentrDistr->GetNbinsX()){
390 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
391 // }
392
393 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
394 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
395 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
396
397 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
398 if(gRandom->Uniform(1.)<bincont)return kTRUE;
399 return kFALSE;
400 }
401
402 if(centDigits*100.<bincont)return kTRUE;
403 return kFALSE;
404
405}
0c22e2ac 406//---------------------------------------------------------------------------
0602cf99 407void AliRDHFCuts::SetupPID(AliVEvent *event) {
408 // Set the PID response object in the AliAODPidHF
409 // in case of old PID sets the TPC dE/dx BB parameterization
e865f048 410
347d8758 411 if(fPidHF){
412 if(fPidHF->GetPidResponse()==0x0){
413 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
414 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
415 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
416 fPidHF->SetPidResponse(pidResp);
417 }
0602cf99 418 if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
347d8758 419 if(fPidHF->GetOldPid()) {
0602cf99 420
421 Bool_t isMC=kFALSE;
422 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
423 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
424
347d8758 425 // pp, from LHC10d onwards
426 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
427 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
428 // pp, 2011 low energy run
429 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
430 fPidHF->SetppLowEn2011(kTRUE);
431 fPidHF->SetOnePad(kFALSE);
432 }
433 // PbPb LHC10h
434 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
435 // MC
436 if(isMC) fPidHF->SetMC(kTRUE);
437 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
438 fPidHF->SetMClowenpp2011(kTRUE);
439 fPidHF->SetBetheBloch();
440 }else{
441 // check that AliPIDResponse object was properly set in case of using OADB
442 if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
0c22e2ac 443 }
0c22e2ac 444 }
0602cf99 445}
446//---------------------------------------------------------------------------
447Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
448 //
449 // Event selection
450 //
451 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
452
453
454
455 fWhyRejection=0;
456 fEvRejectionBits=0;
457 Bool_t accept=kTRUE;
458
0602cf99 459
460 // check if it's MC
461 Bool_t isMC=kFALSE;
462 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
463 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
464
0c22e2ac 465
0602cf99 466 SetupPID(event);
0c22e2ac 467
468 // trigger class
469 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
470 // don't do for MC and for PbPb 2010 data
471 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
68c15f66 472 if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
473 (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
0c22e2ac 474 fWhyRejection=5;
475 fEvRejectionBits+=1<<kNotSelTrigger;
476 accept=kFALSE;
477 }
478 }
479
e865f048 480 // TEMPORARY FIX FOR GetEvent
4640c275 481 Int_t nTracks=((AliAODEvent*)event)->GetNumberOfTracks();
e865f048 482 for(Int_t itr=0; itr<nTracks; itr++){
483 AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
484 tr->SetAODEvent((AliAODEvent*)event);
485 }
486
0c22e2ac 487 // TEMPORARY FIX FOR REFERENCES
488 // Fix references to daughter tracks
489 // if(fFixRefs) {
490 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
491 // fixer->FixReferences((AliAODEvent*)event);
492 // delete fixer;
493 // }
494 //
495
496
497 // physics selection requirements
498 if(fUsePhysicsSelection){
499 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
500 if(!isSelected) {
501 if(accept) fWhyRejection=7;
502 fEvRejectionBits+=1<<kPhysicsSelection;
503 accept=kFALSE;
504 }else{
505 if(fUseOnlyOneTrigger){
506 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
507 if(accept) fWhyRejection=7;
508 fEvRejectionBits+=1<<kPhysicsSelection;
509 accept=kFALSE;
510 }
511 }
512 }
513 }
514
515 // vertex requirements
516
517 const AliVVertex *vertex = event->GetPrimaryVertex();
518
519 if(!vertex){
520 accept=kFALSE;
521 fEvRejectionBits+=1<<kNoVertex;
522 }else{
523 TString title=vertex->GetTitle();
5aa726b3 524 if(title.Contains("Z") && fMinVtxType>1){
525 accept=kFALSE;
526 fEvRejectionBits+=1<<kNoVertex;
527 }
528 else if(title.Contains("3D") && fMinVtxType>2){
529 accept=kFALSE;
530 fEvRejectionBits+=1<<kNoVertex;
0c22e2ac 531 }
532 if(vertex->GetNContributors()<fMinVtxContr){
533 accept=kFALSE;
534 fEvRejectionBits+=1<<kTooFewVtxContrib;
535 }
536 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
537 fEvRejectionBits+=1<<kZVtxOutFid;
538 if(accept) fWhyRejection=6;
539 accept=kFALSE;
540 }
541 }
542
4346ffd2 543 if(fCutOnzVertexSPD>0){
544 const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
545 if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
546 accept=kFALSE;
547 fEvRejectionBits+=1<<kBadSPDVertex;
548 }else{
549 if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
550 fEvRejectionBits+=1<<kZVtxSPDOutFid;
551 if(accept) fWhyRejection=6;
552 accept=kFALSE;
553 }
89b8ddf3 554 if(fCutOnzVertexSPD==2 && vertex){
555 if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
556 fEvRejectionBits+=1<<kZVtxSPDOutFid;
557 if(accept) fWhyRejection=6;
558 accept=kFALSE;
559 }
560 }
4346ffd2 561 }
562 }
0c22e2ac 563
564 // pile-up rejection
565 if(fOptPileup==kRejectPileupEvent){
566 Int_t cutc=(Int_t)fMinContrPileup;
567 Double_t cutz=(Double_t)fMinDzPileup;
568 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
569 if(accept) fWhyRejection=1;
5d3cf93b 570 fEvRejectionBits+=1<<kPileup;
571 accept=kFALSE;
572 }
573 }
574 else if(fOptPileup==kRejectMVPileupEvent){
575 AliAnalysisUtils utils;
576 Bool_t isPUMV = utils.IsPileUpMV(event);
577 if(isPUMV) {
578 if(accept) fWhyRejection=1;
579 fEvRejectionBits+=1<<kPileup;
0c22e2ac 580 accept=kFALSE;
581 }
582 }
583
584 // centrality selection
585 if (fUseCentrality!=kCentOff) {
586 Int_t rejection=IsEventSelectedInCentrality(event);
30399c33 587 Bool_t okCent=kFALSE;
588 if(rejection==0) okCent=kTRUE;
589 if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
590 if(!okCent){
0c22e2ac 591 if(accept) fWhyRejection=rejection;
db0ff379 592 if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
593 else fEvRejectionBits+=1<<kOutsideCentrality;
0c22e2ac 594 accept=kFALSE;
595 }
db0ff379 596
0c22e2ac 597 }
598
8da55fd2 599 // PbPb2011 outliers in tracklets vs. VZERO and centTRK vs. centV0
f63277ef 600 if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
851d1cae 601 if(fRemoveTrackletOutliers){
602 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
603 Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
604 Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
9bf6526e 605 if(ntracklets<1000. && v0cent<cutval){
851d1cae 606 if(accept) fWhyRejection=2;
607 fEvRejectionBits+=1<<kOutsideCentrality;
608 accept=kFALSE;
609 }
610 }
8da55fd2 611 if(fMaxDiffTRKV0Centr>0.){
612 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
613 Double_t trkcent=GetCentrality((AliAODEvent*)event,kCentTRK);
614 if(TMath::Abs(trkcent-v0cent)>fMaxDiffTRKV0Centr){
615 if(accept) fWhyRejection=1;
616 fEvRejectionBits+=1<<kBadTrackV0Correl;
617 accept=kFALSE;
618 }
619 }
851d1cae 620 }
621
180a04cb 622 // Correcting PP2012 flag to remoce tracks crossing SPD misaligned staves for periods 12def
623 if(fApplySPDMisalignedPP2012 && !(event->GetRunNumber()>=195681 && event->GetRunNumber()<=197388)) fApplySPDMisalignedPP2012=false;
624
0c22e2ac 625 return accept;
626}
627//---------------------------------------------------------------------------
09bfd399 628Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const{
0c22e2ac 629 //
630 // Daughter tracks selection
631 //
632 if(!fTrackCuts) return kTRUE;
633
634 Int_t ndaughters = d->GetNDaughters();
635 AliAODVertex *vAOD = d->GetPrimaryVtx();
636 Double_t pos[3],cov[6];
637 vAOD->GetXYZ(pos);
638 vAOD->GetCovarianceMatrix(cov);
639 const AliESDVertex vESD(pos,cov,100.,100);
640
641 Bool_t retval=kTRUE;
642
643 for(Int_t idg=0; idg<ndaughters; idg++) {
644 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
645 if(!dgTrack) {retval = kFALSE; continue;}
646 //printf("charge %d\n",dgTrack->Charge());
647 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
648
649 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
650 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
651
652 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
653 }
654
655 return retval;
656}
657//---------------------------------------------------------------------------
09bfd399 658Bool_t AliRDHFCuts::CheckPtDepCrossedRows(TString rows,Bool_t print) const {
659 //
660 // Check the correctness of the string syntax
661 //
662 Bool_t retval=kTRUE;
663
664 if(!rows.Contains("pt")) {
665 if(print) AliError("string must contain \"pt\"");
666 retval= kFALSE;
667 }
668 return retval;
669}
670//---------------------------------------------------------------------------
671void AliRDHFCuts::SetMinCrossedRowsTPCPtDep(const char *rows){
672 //
673 //Create the TFormula from TString for TPC crossed rows pT dependent cut
674 //
675
676
677 // setting data member that describes the TPC crossed rows pT dependent cut
678 fCutMinCrossedRowsTPCPtDep = rows;
679
680 // creating TFormula from TString
681 if(f1CutMinNCrossedRowsTPCPtDep){
682 delete f1CutMinNCrossedRowsTPCPtDep;
683 // resetting TFormula
684 f1CutMinNCrossedRowsTPCPtDep = 0;
685 }
686 if(!CheckPtDepCrossedRows(rows,kTRUE))return;
687
688 TString tmp(rows);
689 tmp.ReplaceAll("pt","x");
690 f1CutMinNCrossedRowsTPCPtDep = new TFormula("f1CutMinNCrossedRowsTPCPtDep",tmp.Data());
691
692
693}
694//---------------------------------------------------------------------------
695Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const{
0c22e2ac 696 //
697 // Convert to ESDtrack, relate to vertex and check cuts
698 //
699 if(!cuts) return kTRUE;
700
701 if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
702
0c22e2ac 703
704 // convert to ESD track here
705 AliESDtrack esdTrack(track);
706 // set the TPC cluster info
707 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
708 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
709 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
710 // needed to calculate the impact parameters
09bfd399 711 esdTrack.RelateToVertex(primary,0.,3.);
fdd5a897 712
09bfd399 713 //applying ESDtrackCut
714 if(!cuts->IsSelected(&esdTrack)) return kFALSE;
fdd5a897 715
09bfd399 716 //appliyng kink rejection
fdd5a897 717 if(fKinkReject){
718 AliAODVertex *maybeKink=track->GetProdVertex();
09bfd399 719 if(maybeKink->GetType()==AliAODVertex::kKink) return kFALSE;
720 }
721
722 //appliyng TPC crossed rows pT dependent cut
723 if(f1CutMinNCrossedRowsTPCPtDep){
724 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
725 if(nCrossedRowsTPC<f1CutMinNCrossedRowsTPCPtDep->Eval(esdTrack.Pt())) return kFALSE;
fdd5a897 726 }
09bfd399 727
728 //appliyng NTPCcls/NTPCcrossedRows cut
729 if(fCutRatioClsOverCrossRowsTPC){
730 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
731 Float_t nClustersTPC = esdTrack.GetTPCNcls();
732 if(nCrossedRowsTPC!=0){
733 Float_t ratio = nClustersTPC/nCrossedRowsTPC;
734 if(ratio<fCutRatioClsOverCrossRowsTPC) return kFALSE;
735 }
736 else return kFALSE;
737 }
738
739 //appliyng TPCsignalN/NTPCcrossedRows cut
740 if(fCutRatioSignalNOverCrossRowsTPC){
741 Float_t nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
742 Float_t nTPCsignal = esdTrack.GetTPCsignalN();
743 if(nCrossedRowsTPC!=0){
744 Float_t ratio = nTPCsignal/nCrossedRowsTPC;
745 if(ratio<fCutRatioSignalNOverCrossRowsTPC) return kFALSE;
746 }
747 else return kFALSE;
748 }
749
0c22e2ac 750
751 if(fOptPileup==kRejectTracksFromPileupVertex){
752 // to be implemented
753 // we need either to have here the AOD Event,
754 // or to have the pileup vertex object
755 }
851d1cae 756
757 if(fApplySPDDeadPbPb2011){
758 Bool_t deadSPDLay1PbPb2011[20][4]={
759 {kTRUE,kTRUE,kTRUE,kTRUE},
760 {kTRUE,kTRUE,kTRUE,kTRUE},
761 {kTRUE,kTRUE,kTRUE,kTRUE},
762 {kTRUE,kTRUE,kTRUE,kTRUE},
763 {kTRUE,kTRUE,kTRUE,kTRUE},
764 {kFALSE,kFALSE,kTRUE,kTRUE},
765 {kTRUE,kTRUE,kFALSE,kFALSE},
766 {kTRUE,kTRUE,kTRUE,kTRUE},
767 {kFALSE,kFALSE,kTRUE,kTRUE},
768 {kTRUE,kTRUE,kTRUE,kTRUE},
769 {kTRUE,kTRUE,kFALSE,kFALSE},
770 {kTRUE,kTRUE,kTRUE,kTRUE},
771 {kFALSE,kFALSE,kFALSE,kFALSE},
772 {kFALSE,kFALSE,kTRUE,kTRUE},
773 {kFALSE,kFALSE,kFALSE,kFALSE},
774 {kFALSE,kFALSE,kFALSE,kFALSE},
775 {kTRUE,kTRUE,kTRUE,kTRUE},
776 {kTRUE,kTRUE,kFALSE,kFALSE},
777 {kFALSE,kFALSE,kFALSE,kFALSE},
778 {kFALSE,kFALSE,kFALSE,kFALSE}
779 };
780 Bool_t deadSPDLay2PbPb2011[40][4]={
781 {kTRUE,kTRUE,kTRUE,kTRUE},
782 {kTRUE,kTRUE,kTRUE,kTRUE},
783 {kTRUE,kTRUE,kTRUE,kTRUE},
784 {kTRUE,kTRUE,kTRUE,kTRUE},
785 {kTRUE,kTRUE,kTRUE,kTRUE},
786 {kTRUE,kTRUE,kTRUE,kTRUE},
787 {kTRUE,kTRUE,kTRUE,kTRUE},
788 {kTRUE,kTRUE,kTRUE,kTRUE},
789 {kTRUE,kTRUE,kTRUE,kTRUE},
790 {kTRUE,kTRUE,kTRUE,kTRUE},
791 {kTRUE,kTRUE,kTRUE,kTRUE},
792 {kTRUE,kTRUE,kTRUE,kTRUE},
793 {kFALSE,kFALSE,kFALSE,kFALSE},
794 {kFALSE,kFALSE,kTRUE,kTRUE},
795 {kTRUE,kTRUE,kTRUE,kTRUE},
796 {kTRUE,kTRUE,kTRUE,kTRUE},
797 {kTRUE,kTRUE,kFALSE,kFALSE},
798 {kTRUE,kTRUE,kTRUE,kTRUE},
799 {kTRUE,kTRUE,kTRUE,kTRUE},
800 {kTRUE,kTRUE,kTRUE,kTRUE},
801 {kFALSE,kFALSE,kFALSE,kFALSE},
802 {kFALSE,kFALSE,kFALSE,kFALSE},
803 {kTRUE,kTRUE,kTRUE,kTRUE},
804 {kTRUE,kTRUE,kTRUE,kTRUE},
805 {kFALSE,kFALSE,kFALSE,kFALSE},
806 {kFALSE,kFALSE,kFALSE,kFALSE},
807 {kTRUE,kTRUE,kTRUE,kTRUE},
808 {kTRUE,kTRUE,kTRUE,kTRUE},
809 {kFALSE,kFALSE,kFALSE,kFALSE},
810 {kFALSE,kFALSE,kFALSE,kFALSE},
811 {kFALSE,kFALSE,kFALSE,kFALSE},
812 {kFALSE,kFALSE,kFALSE,kFALSE},
813 {kTRUE,kTRUE,kTRUE,kTRUE},
814 {kTRUE,kTRUE,kTRUE,kTRUE},
815 {kTRUE,kTRUE,kFALSE,kFALSE},
816 {kTRUE,kTRUE,kTRUE,kTRUE},
817 {kFALSE,kFALSE,kFALSE,kFALSE},
818 {kFALSE,kFALSE,kFALSE,kFALSE},
819 {kFALSE,kFALSE,kFALSE,kFALSE},
820 {kFALSE,kFALSE,kFALSE,kFALSE}
821 };
822 Double_t xyz1[3],xyz2[3];
823 esdTrack.GetXYZAt(3.9,0.,xyz1);
824 esdTrack.GetXYZAt(7.6,0.,xyz2);
825 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
826 if(phi1<0) phi1+=2*TMath::Pi();
a0dabf3d 827 Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
851d1cae 828 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
829 if(phi2<0) phi2+=2*TMath::Pi();
a0dabf3d 830 Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
831 Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
832 Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
851d1cae 833 Bool_t lay1ok=kFALSE;
834 if(mod1>=0 && mod1<4 && lad1<20){
835 lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
836 }
837 Bool_t lay2ok=kFALSE;
838 if(mod2>=0 && mod2<4 && lad2<40){
839 lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
840 }
09bfd399 841 if(!lay1ok && !lay2ok) return kFALSE;
851d1cae 842 }
843
180a04cb 844 if(fApplySPDMisalignedPP2012) {
845 // Cut tracks crossing the SPD at 5.6<phi<2pi
846 Double_t xyz1[3],xyz2[3];
847 esdTrack.GetXYZAt(3.9,0.,xyz1);
848 esdTrack.GetXYZAt(7.6,0.,xyz2);
849 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
850 if(phi1<0) phi1+=2*TMath::Pi();
851 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
852 if(phi2<0) phi2+=2*TMath::Pi();
853 Bool_t lay1ok=kTRUE;
854 if(phi1>5.6 && phi1<2.*TMath::Pi()) lay1ok=kFALSE;
855 Bool_t lay2ok=kTRUE;
856 if(phi2>5.6 && phi2<2.*TMath::Pi()) lay2ok=kFALSE;
857 if(!lay1ok || !lay2ok) return kFALSE;
858 }
859
09bfd399 860 return kTRUE;
0c22e2ac 861}
862//---------------------------------------------------------------------------
863void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
864 // Set the pt bins
865
866 if(fPtBinLimits) {
867 delete [] fPtBinLimits;
868 fPtBinLimits = NULL;
869 printf("Changing the pt bins\n");
870 }
871
872 if(nPtBinLimits != fnPtBins+1){
873 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
874 SetNPtBins(nPtBinLimits-1);
875 }
876
877 fnPtBinLimits = nPtBinLimits;
878 SetGlobalIndex();
879 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
880 fPtBinLimits = new Float_t[fnPtBinLimits];
881 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
882
883 return;
884}
885//---------------------------------------------------------------------------
886void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
887 // Set the variable names
888
889 if(fVarNames) {
890 delete [] fVarNames;
891 fVarNames = NULL;
892 //printf("Changing the variable names\n");
893 }
894 if(nVars!=fnVars){
895 printf("Wrong number of variables: it has to be %d\n",fnVars);
896 return;
897 }
898 //fnVars=nVars;
899 fVarNames = new TString[nVars];
900 fIsUpperCut = new Bool_t[nVars];
901 for(Int_t iv=0; iv<nVars; iv++) {
902 fVarNames[iv] = varNames[iv];
903 fIsUpperCut[iv] = isUpperCut[iv];
904 }
905
906 return;
907}
908//---------------------------------------------------------------------------
909void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
910 // Set the variables to be used for cuts optimization
911
912 if(fVarsForOpt) {
913 delete [] fVarsForOpt;
914 fVarsForOpt = NULL;
915 //printf("Changing the variables for cut optimization\n");
916 }
917
918 if(nVars==0){//!=fnVars) {
919 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
920 return;
921 }
922
923 fnVarsForOpt = 0;
924 fVarsForOpt = new Bool_t[fnVars];
925 for(Int_t iv=0; iv<fnVars; iv++) {
926 fVarsForOpt[iv]=forOpt[iv];
927 if(fVarsForOpt[iv]) fnVarsForOpt++;
928 }
929
930 return;
931}
932
933//---------------------------------------------------------------------------
934void AliRDHFCuts::SetUseCentrality(Int_t flag) {
935 //
936 // set centrality estimator
937 //
938 fUseCentrality=flag;
939 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
940
941 return;
942}
943
944
945//---------------------------------------------------------------------------
946void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
947 //
948 // store the cuts
949 //
950 if(nVars!=fnVars) {
951 printf("Wrong number of variables: it has to be %d\n",fnVars);
952 AliFatal("exiting");
953 }
954 if(nPtBins!=fnPtBins) {
955 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
956 AliFatal("exiting");
957 }
958
959 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
960
961
962 for(Int_t iv=0; iv<fnVars; iv++) {
963
964 for(Int_t ib=0; ib<fnPtBins; ib++) {
965
966 //check
967 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
968 cout<<"Overflow, exit..."<<endl;
969 return;
970 }
971
972 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
973
974 }
975 }
976 return;
977}
978//---------------------------------------------------------------------------
979void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
980 //
981 // store the cuts
982 //
983 if(glIndex != fGlobalIndex){
984 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
985 AliFatal("exiting");
986 }
987 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
988
989 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
990 fCutsRD[iGl] = cutsRDGlob[iGl];
991 }
992 return;
993}
994//---------------------------------------------------------------------------
995void AliRDHFCuts::PrintAll() const {
996 //
997 // print all cuts values
998 //
999
1000 printf("Minimum vtx type %d\n",fMinVtxType);
1001 printf("Minimum vtx contr %d\n",fMinVtxContr);
1002 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
1003 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
e865f048 1004 printf("Use PID %d OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
0c22e2ac 1005 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
1006 printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
1007 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
1008 if(fOptPileup==1) printf(" -- Reject pileup event");
1009 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
1010 if(fUseCentrality>0) {
1011 TString estimator="";
1012 if(fUseCentrality==1) estimator = "V0";
1013 if(fUseCentrality==2) estimator = "Tracks";
1014 if(fUseCentrality==3) estimator = "Tracklets";
1015 if(fUseCentrality==4) estimator = "SPD clusters outer";
cc2b144d 1016 if(fUseCentrality==5) estimator = "ZNA";
1017 if(fUseCentrality==6) estimator = "ZPA";
1018 if(fUseCentrality==7) estimator = "V0A";
09bfd399 1019 printf("Centrality class considered: %.1f-%.1f, estimated with %s\n",fMinCentrality,fMaxCentrality,estimator.Data());
0c22e2ac 1020 }
1021 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
1022
09bfd399 1023 if(fCutRatioClsOverCrossRowsTPC) printf("N TPC Clusters > %f N TPC Crossed Rows\n", fCutRatioClsOverCrossRowsTPC);
1024 if(fCutRatioSignalNOverCrossRowsTPC) printf("N TPC Points for dE/dx > %f N TPC Crossed Rows\n", fCutRatioSignalNOverCrossRowsTPC);
1025 if(f1CutMinNCrossedRowsTPCPtDep) printf("N TPC Crossed Rows pT-dependent cut: %s\n", fCutMinCrossedRowsTPCPtDep.Data());
1026
0c22e2ac 1027 if(fVarNames){
1028 cout<<"Array of variables"<<endl;
1029 for(Int_t iv=0;iv<fnVars;iv++){
1030 cout<<fVarNames[iv]<<"\t";
1031 }
1032 cout<<endl;
1033 }
1034 if(fVarsForOpt){
1035 cout<<"Array of optimization"<<endl;
1036 for(Int_t iv=0;iv<fnVars;iv++){
1037 cout<<fVarsForOpt[iv]<<"\t";
1038 }
1039 cout<<endl;
1040 }
1041 if(fIsUpperCut){
1042 cout<<"Array of upper/lower cut"<<endl;
1043 for(Int_t iv=0;iv<fnVars;iv++){
1044 cout<<fIsUpperCut[iv]<<"\t";
1045 }
1046 cout<<endl;
1047 }
1048 if(fPtBinLimits){
1049 cout<<"Array of ptbin limits"<<endl;
1050 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
1051 cout<<fPtBinLimits[ib]<<"\t";
1052 }
1053 cout<<endl;
1054 }
1055 if(fCutsRD){
1056 cout<<"Matrix of cuts"<<endl;
1057 for(Int_t iv=0;iv<fnVars;iv++){
1058 for(Int_t ib=0;ib<fnPtBins;ib++){
1059 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
1060 }
1061 cout<<endl;
1062 }
1063 cout<<endl;
1064 }
5de70ece 1065 if(fPidHF) fPidHF->PrintAll();
0c22e2ac 1066 return;
1067}
1068
1069//--------------------------------------------------------------------------
1070void AliRDHFCuts::PrintTrigger() const{
b8e9fa32 1071 // print the trigger selection
1072
1073 printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
0c22e2ac 1074
1075 cout<<" Trigger selection pattern: ";
1076 if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
1077 if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
1078 if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
1079 if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
1080 if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
1081 if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
1082 if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
1083 if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
1084 if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
1085 if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
1086 cout << endl<< endl;
1087
1088}
1089
1090//---------------------------------------------------------------------------
1091void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
1092 //
1093 // get the cuts
1094 //
1095
1096 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
1097
1098
1099 Int_t iv,ib;
1100 if(!cutsRD) {
1101 //cout<<"Initialization..."<<endl;
1102 cutsRD=new Float_t*[fnVars];
1103 for(iv=0; iv<fnVars; iv++) {
1104 cutsRD[iv] = new Float_t[fnPtBins];
1105 }
1106 }
1107
1108 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
1109 GetVarPtIndex(iGlobal,iv,ib);
1110 cutsRD[iv][ib] = fCutsRD[iGlobal];
1111 }
1112
1113 return;
1114}
1115
1116//---------------------------------------------------------------------------
1117Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
1118 //
1119 // give the global index from variable and pt bin
1120 //
1121 return iPtBin*fnVars+iVar;
1122}
1123
1124//---------------------------------------------------------------------------
1125void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
1126 //
1127 //give the index of the variable and of the pt bin from the global index
1128 //
1129 iPtBin=(Int_t)iGlob/fnVars;
1130 iVar=iGlob%fnVars;
1131
1132 return;
1133}
1134
1135//---------------------------------------------------------------------------
1136Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1137 //
1138 //give the pt bin where the pt lies.
1139 //
1140 Int_t ptbin=-1;
1141 if(pt<fPtBinLimits[0])return ptbin;
1142 for (Int_t i=0;i<fnPtBins;i++){
1143 if(pt<fPtBinLimits[i+1]) {
1144 ptbin=i;
1145 break;
1146 }
1147 }
1148 return ptbin;
1149}
1150//-------------------------------------------------------------------
1151Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1152 //
1153 // Give the value of cut set for the variable iVar and the pt bin iPtBin
1154 //
1155 if(!fCutsRD){
1156 cout<<"Cuts not iniziaisez yet"<<endl;
1157 return 0;
1158 }
1159 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1160}
1161//-------------------------------------------------------------------
1162Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1163 //
1164 // Get centrality percentile
1165 //
1166
1167 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1168 if(mcArray) {fUseAOD049=kFALSE;}
1169
0a918d8d 1170 AliAODHeader *header=dynamic_cast<AliAODHeader*>(aodEvent->GetHeader());
1171 if(!header) AliFatal("Not a standard AOD");
0c22e2ac 1172 AliCentrality *centrality=header->GetCentralityP();
1173 Float_t cent=-999.;
1174 Bool_t isSelRun=kFALSE;
1175 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1176 if(!centrality) return cent;
1177 else{
1178 if (estimator==kCentV0M){
1179 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1180 if(cent<0){
1181 Int_t quality = centrality->GetQuality();
4346ffd2 1182 if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
0c22e2ac 1183 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1184 }else{
1185 Int_t runnum=aodEvent->GetRunNumber();
1186 for(Int_t ir=0;ir<5;ir++){
1187 if(runnum==selRun[ir]){
1188 isSelRun=kTRUE;
1189 break;
1190 }
1191 }
1192 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1193 }
1194 }
1195
1196 //temporary fix for AOD049 outliers
1197 if(fUseAOD049&&cent>=0){
1198 Float_t v0=0;
1199 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1200 v0+=aodV0->GetMTotV0A();
1201 v0+=aodV0->GetMTotV0C();
1202 if(cent==0&&v0<19500)return -1;//filtering issue
1203 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1204 Float_t val= 1.30552 + 0.147931 * v0;
1205 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};
1206 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1207 }
1208 }
1209 else {
cc2b144d 1210 if (estimator==kCentTRK) {
0c22e2ac 1211 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1212 if(cent<0){
1213 Int_t quality = centrality->GetQuality();
1214 if(quality<=1){
1215 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1216 }else{
1217 Int_t runnum=aodEvent->GetRunNumber();
1218 for(Int_t ir=0;ir<5;ir++){
1219 if(runnum==selRun[ir]){
1220 isSelRun=kTRUE;
1221 break;
1222 }
1223 }
1224 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1225 }
1226 }
cc2b144d 1227 }
0c22e2ac 1228 else{
1229 if (estimator==kCentTKL){
1230 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1231 if(cent<0){
1232 Int_t quality = centrality->GetQuality();
1233 if(quality<=1){
1234 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1235 }else{
1236 Int_t runnum=aodEvent->GetRunNumber();
1237 for(Int_t ir=0;ir<5;ir++){
1238 if(runnum==selRun[ir]){
1239 isSelRun=kTRUE;
1240 break;
1241 }
1242 }
1243 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1244 }
1245 }
1246 }
1247 else{
1248 if (estimator==kCentCL1){
1249 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1250 if(cent<0){
1251 Int_t quality = centrality->GetQuality();
1252 if(quality<=1){
1253 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1254 }else{
1255 Int_t runnum=aodEvent->GetRunNumber();
1256 for(Int_t ir=0;ir<5;ir++){
1257 if(runnum==selRun[ir]){
1258 isSelRun=kTRUE;
1259 break;
1260 }
1261 }
1262 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1263 }
1264 }
1265 }
cc2b144d 1266 else{
1267 if (estimator==kCentZNA){
1268 cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1269 if(cent<0){
1270 Int_t quality = centrality->GetQuality();
1271 if(quality<=1){
1272 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1273 }else{
1274 Int_t runnum=aodEvent->GetRunNumber();
1275 for(Int_t ir=0;ir<5;ir++){
1276 if(runnum==selRun[ir]){
1277 isSelRun=kTRUE;
1278 break;
1279 }
1280 }
1281 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1282 }
1283 }
1284 }
1285 else{
1286 if (estimator==kCentZPA){
1287 cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1288 if(cent<0){
1289 Int_t quality = centrality->GetQuality();
1290 if(quality<=1){
1291 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1292 }else{
1293 Int_t runnum=aodEvent->GetRunNumber();
1294 for(Int_t ir=0;ir<5;ir++){
1295 if(runnum==selRun[ir]){
1296 isSelRun=kTRUE;
1297 break;
1298 }
1299 }
1300 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1301 }
1302 }
1303 }
1304 else{
1305 if (estimator==kCentV0A){
1306 cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1307 if(cent<0){
1308 Int_t quality = centrality->GetQuality();
1309 if(quality<=1){
1310 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1311 }else{
1312 Int_t runnum=aodEvent->GetRunNumber();
1313 for(Int_t ir=0;ir<5;ir++){
1314 if(runnum==selRun[ir]){
1315 isSelRun=kTRUE;
1316 break;
1317 }
1318 }
1319 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1320 }
1321 }
1322 }
0c22e2ac 1323 else {
1324 AliWarning("Centrality estimator not valid");
1325
1326 }
1327 }
cc2b144d 1328 }
1329 }
1330 }
1331 }
1332 }
0c22e2ac 1333 }
1334 return cent;
1335}
1336//-------------------------------------------------------------------
1337Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1338 //
1339 // Compare two cuts objects
1340 //
1341
1342 Bool_t areEqual=kTRUE;
1343
1344 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1345
1346 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1347
1348 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1349
1350 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1351
1352 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1353
1354 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1355 if(fTrackCuts){
1356 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1357
1358 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1359
1360 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1361
1362 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;}
1363 }
1364
1365 if(fCutsRD) {
1366 for(Int_t iv=0;iv<fnVars;iv++) {
1367 for(Int_t ib=0;ib<fnPtBins;ib++) {
1368 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1369 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1370 areEqual=kFALSE;
1371 }
1372 }
1373 }
1374 }
1375
1376 return areEqual;
1377}
1378//---------------------------------------------------------------------------
1379void AliRDHFCuts::MakeTable() const {
1380 //
1381 // print cuts values in table format
1382 //
1383
1384 TString ptString = "pT range";
1385 if(fVarNames && fPtBinLimits && fCutsRD){
1386 TString firstLine(Form("* %-15s",ptString.Data()));
1387 for (Int_t ivar=0; ivar<fnVars; ivar++){
1388 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1389 if (ivar == fnVars){
1390 firstLine+="*\n";
1391 }
1392 }
1393 Printf("%s",firstLine.Data());
1394
1395 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1396 TString line;
1397 if (ipt==fnPtBins-1){
1398 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1399 }
1400 else{
1401 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1402 }
1403 for (Int_t ivar=0; ivar<fnVars; ivar++){
1404 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1405 }
1406 Printf("%s",line.Data());
1407 }
1408
1409 }
1410
1411
1412 return;
1413}
1414//--------------------------------------------------------------------------
1415Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1416 AliAODEvent *aod) const
1417{
1418 //
1419 // Recalculate primary vertex without daughters
1420 //
1421
1422 if(!aod) {
1423 AliError("Can not remove daughters from vertex without AOD event");
1424 return 0;
1425 }
1426
1427 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1428 if(!recvtx){
1429 AliDebug(2,"Removal of daughter tracks failed");
1430 return kFALSE;
1431 }
1432
1433
1434 //set recalculed primary vertex
1435 d->SetOwnPrimaryVtx(recvtx);
1436 delete recvtx;
1437
1438 return kTRUE;
1439}
1440//--------------------------------------------------------------------------
1441Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1442{
1443 //
1444 // Recalculate primary vertex without daughters
1445 //
1446
1447 if(!aod) {
1448 AliError("Can not get MC vertex without AOD event");
1449 return kFALSE;
1450 }
1451
1452 // load MC header
1453 AliAODMCHeader *mcHeader =
1454 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1455 if(!mcHeader) {
1456 AliError("Can not get MC vertex without AODMCHeader event");
1457 return kFALSE;
1458 }
1459 Double_t pos[3];
1460 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1461 mcHeader->GetVertex(pos);
1462 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1463
1464 if(!recvtx){
1465 AliDebug(2,"Removal of daughter tracks failed");
1466 return kFALSE;
1467 }
1468
1469 //set recalculed primary vertex
1470 d->SetOwnPrimaryVtx(recvtx);
1471
1472 d->RecalculateImpPars(recvtx,aod);
1473
1474 delete recvtx;
1475
1476 return kTRUE;
1477}
1478//--------------------------------------------------------------------------
1479void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1480 AliAODEvent *aod,
1481 AliAODVertex *origownvtx) const
1482{
1483 //
1484 // Clean-up own primary vertex if needed
1485 //
1486
1487 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1488 d->UnsetOwnPrimaryVtx();
1489 if(origownvtx) {
1490 d->SetOwnPrimaryVtx(origownvtx);
1491 delete origownvtx; origownvtx=NULL;
1492 }
1493 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1494 } else {
1495 if(origownvtx) {
1496 delete origownvtx; origownvtx=NULL;
1497 }
1498 }
1499 return;
1500}
1501//--------------------------------------------------------------------------
1502Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1503{
1504 //
1505 // Checks if this candidate is matched to MC signal
1506 //
1507
1508 if(!aod) return kFALSE;
1509
1510 // get the MC array
1511 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1512
1513 if(!mcArray) return kFALSE;
1514
1515 // try to match
1516 Int_t label = d->MatchToMC(pdg,mcArray);
1517
1518 if(label>=0) {
1519 //printf("MATCH!\n");
1520 return kTRUE;
1521 }
1522
1523 return kFALSE;
1524}
1525
1526
4bdc0315 1527//--------------------------------------------------------------------------
1528Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1529 // recompute event primary vertex from AOD tracks
1530
1531 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1532 vertexer->SetITSMode();
1533 vertexer->SetMinClusters(3);
1534
1535 AliAODVertex* pvtx=event->GetPrimaryVertex();
1536 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1537 Float_t diamondcovxy[3];
1538 event->GetDiamondCovXY(diamondcovxy);
1539 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1540 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1541 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1542 vertexer->SetVtxStart(diamond);
1543 delete diamond; diamond=NULL;
1544 }
1545
1546 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1547 if(!vertexESD) return kFALSE;
1548 if(vertexESD->GetNContributors()<=0) {
1549 //AliDebug(2,"vertexing failed");
1550 delete vertexESD; vertexESD=NULL;
1551 return kFALSE;
1552 }
1553 delete vertexer; vertexer=NULL;
1554
1555 // convert to AliAODVertex
1556 Double_t pos[3],cov[6],chi2perNDF;
1557 vertexESD->GetXYZ(pos); // position
1558 vertexESD->GetCovMatrix(cov); //covariance matrix
1559 chi2perNDF = vertexESD->GetChi2toNDF();
1560 delete vertexESD; vertexESD=NULL;
1561
1562 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1563 pvtx->SetChi2perNDF(chi2perNDF);
1564 pvtx->SetCovMatrix(cov);
1565
1566 return kTRUE;
1567}