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