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