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