]> 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
consolidate zero-length arrays (aka struct hack)
[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
491 Int_t nTracks=((AliAODEvent*)event)->GetNTracks();
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
1181 AliAODHeader *header=aodEvent->GetHeader();
1182 AliCentrality *centrality=header->GetCentralityP();
1183 Float_t cent=-999.;
1184 Bool_t isSelRun=kFALSE;
1185 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1186 if(!centrality) return cent;
1187 else{
1188 if (estimator==kCentV0M){
1189 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1190 if(cent<0){
1191 Int_t quality = centrality->GetQuality();
4346ffd2 1192 if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
0c22e2ac 1193 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1194 }else{
1195 Int_t runnum=aodEvent->GetRunNumber();
1196 for(Int_t ir=0;ir<5;ir++){
1197 if(runnum==selRun[ir]){
1198 isSelRun=kTRUE;
1199 break;
1200 }
1201 }
1202 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1203 }
1204 }
1205
1206 //temporary fix for AOD049 outliers
1207 if(fUseAOD049&&cent>=0){
1208 Float_t v0=0;
1209 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1210 v0+=aodV0->GetMTotV0A();
1211 v0+=aodV0->GetMTotV0C();
1212 if(cent==0&&v0<19500)return -1;//filtering issue
1213 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1214 Float_t val= 1.30552 + 0.147931 * v0;
1215 Float_t tklSigma[101]={176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86, 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654, 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334, 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224, 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255, 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398, 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235, 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504, 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544};
1216 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1217 }
1218 }
1219 else {
cc2b144d 1220 if (estimator==kCentTRK) {
0c22e2ac 1221 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1222 if(cent<0){
1223 Int_t quality = centrality->GetQuality();
1224 if(quality<=1){
1225 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1226 }else{
1227 Int_t runnum=aodEvent->GetRunNumber();
1228 for(Int_t ir=0;ir<5;ir++){
1229 if(runnum==selRun[ir]){
1230 isSelRun=kTRUE;
1231 break;
1232 }
1233 }
1234 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1235 }
1236 }
cc2b144d 1237 }
0c22e2ac 1238 else{
1239 if (estimator==kCentTKL){
1240 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1241 if(cent<0){
1242 Int_t quality = centrality->GetQuality();
1243 if(quality<=1){
1244 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1245 }else{
1246 Int_t runnum=aodEvent->GetRunNumber();
1247 for(Int_t ir=0;ir<5;ir++){
1248 if(runnum==selRun[ir]){
1249 isSelRun=kTRUE;
1250 break;
1251 }
1252 }
1253 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1254 }
1255 }
1256 }
1257 else{
1258 if (estimator==kCentCL1){
1259 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1260 if(cent<0){
1261 Int_t quality = centrality->GetQuality();
1262 if(quality<=1){
1263 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1264 }else{
1265 Int_t runnum=aodEvent->GetRunNumber();
1266 for(Int_t ir=0;ir<5;ir++){
1267 if(runnum==selRun[ir]){
1268 isSelRun=kTRUE;
1269 break;
1270 }
1271 }
1272 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1273 }
1274 }
1275 }
cc2b144d 1276 else{
1277 if (estimator==kCentZNA){
1278 cent=(Float_t)(centrality->GetCentralityPercentile("ZNA"));
1279 if(cent<0){
1280 Int_t quality = centrality->GetQuality();
1281 if(quality<=1){
1282 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1283 }else{
1284 Int_t runnum=aodEvent->GetRunNumber();
1285 for(Int_t ir=0;ir<5;ir++){
1286 if(runnum==selRun[ir]){
1287 isSelRun=kTRUE;
1288 break;
1289 }
1290 }
1291 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZNA");
1292 }
1293 }
1294 }
1295 else{
1296 if (estimator==kCentZPA){
1297 cent=(Float_t)(centrality->GetCentralityPercentile("ZPA"));
1298 if(cent<0){
1299 Int_t quality = centrality->GetQuality();
1300 if(quality<=1){
1301 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1302 }else{
1303 Int_t runnum=aodEvent->GetRunNumber();
1304 for(Int_t ir=0;ir<5;ir++){
1305 if(runnum==selRun[ir]){
1306 isSelRun=kTRUE;
1307 break;
1308 }
1309 }
1310 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("ZPA");
1311 }
1312 }
1313 }
1314 else{
1315 if (estimator==kCentV0A){
1316 cent=(Float_t)(centrality->GetCentralityPercentile("V0A"));
1317 if(cent<0){
1318 Int_t quality = centrality->GetQuality();
1319 if(quality<=1){
1320 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1321 }else{
1322 Int_t runnum=aodEvent->GetRunNumber();
1323 for(Int_t ir=0;ir<5;ir++){
1324 if(runnum==selRun[ir]){
1325 isSelRun=kTRUE;
1326 break;
1327 }
1328 }
1329 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0A");
1330 }
1331 }
1332 }
0c22e2ac 1333 else {
1334 AliWarning("Centrality estimator not valid");
1335
1336 }
1337 }
cc2b144d 1338 }
1339 }
1340 }
1341 }
1342 }
0c22e2ac 1343 }
1344 return cent;
1345}
1346//-------------------------------------------------------------------
1347Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1348 //
1349 // Compare two cuts objects
1350 //
1351
1352 Bool_t areEqual=kTRUE;
1353
1354 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1355
1356 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1357
1358 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1359
1360 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1361
1362 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1363
1364 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1365 if(fTrackCuts){
1366 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1367
1368 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1369
1370 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1371
1372 if(fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)!=obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)) {printf("ClusterReq SPD %d %d\n",fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD),obj->fTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD)); areEqual=kFALSE;}
1373 }
1374
1375 if(fCutsRD) {
1376 for(Int_t iv=0;iv<fnVars;iv++) {
1377 for(Int_t ib=0;ib<fnPtBins;ib++) {
1378 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1379 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1380 areEqual=kFALSE;
1381 }
1382 }
1383 }
1384 }
1385
1386 return areEqual;
1387}
1388//---------------------------------------------------------------------------
1389void AliRDHFCuts::MakeTable() const {
1390 //
1391 // print cuts values in table format
1392 //
1393
1394 TString ptString = "pT range";
1395 if(fVarNames && fPtBinLimits && fCutsRD){
1396 TString firstLine(Form("* %-15s",ptString.Data()));
1397 for (Int_t ivar=0; ivar<fnVars; ivar++){
1398 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1399 if (ivar == fnVars){
1400 firstLine+="*\n";
1401 }
1402 }
1403 Printf("%s",firstLine.Data());
1404
1405 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1406 TString line;
1407 if (ipt==fnPtBins-1){
1408 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1409 }
1410 else{
1411 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1412 }
1413 for (Int_t ivar=0; ivar<fnVars; ivar++){
1414 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1415 }
1416 Printf("%s",line.Data());
1417 }
1418
1419 }
1420
1421
1422 return;
1423}
1424//--------------------------------------------------------------------------
1425Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1426 AliAODEvent *aod) const
1427{
1428 //
1429 // Recalculate primary vertex without daughters
1430 //
1431
1432 if(!aod) {
1433 AliError("Can not remove daughters from vertex without AOD event");
1434 return 0;
1435 }
1436
1437 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1438 if(!recvtx){
1439 AliDebug(2,"Removal of daughter tracks failed");
1440 return kFALSE;
1441 }
1442
1443
1444 //set recalculed primary vertex
1445 d->SetOwnPrimaryVtx(recvtx);
1446 delete recvtx;
1447
1448 return kTRUE;
1449}
1450//--------------------------------------------------------------------------
1451Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1452{
1453 //
1454 // Recalculate primary vertex without daughters
1455 //
1456
1457 if(!aod) {
1458 AliError("Can not get MC vertex without AOD event");
1459 return kFALSE;
1460 }
1461
1462 // load MC header
1463 AliAODMCHeader *mcHeader =
1464 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1465 if(!mcHeader) {
1466 AliError("Can not get MC vertex without AODMCHeader event");
1467 return kFALSE;
1468 }
1469 Double_t pos[3];
1470 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1471 mcHeader->GetVertex(pos);
1472 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1473
1474 if(!recvtx){
1475 AliDebug(2,"Removal of daughter tracks failed");
1476 return kFALSE;
1477 }
1478
1479 //set recalculed primary vertex
1480 d->SetOwnPrimaryVtx(recvtx);
1481
1482 d->RecalculateImpPars(recvtx,aod);
1483
1484 delete recvtx;
1485
1486 return kTRUE;
1487}
1488//--------------------------------------------------------------------------
1489void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1490 AliAODEvent *aod,
1491 AliAODVertex *origownvtx) const
1492{
1493 //
1494 // Clean-up own primary vertex if needed
1495 //
1496
1497 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1498 d->UnsetOwnPrimaryVtx();
1499 if(origownvtx) {
1500 d->SetOwnPrimaryVtx(origownvtx);
1501 delete origownvtx; origownvtx=NULL;
1502 }
1503 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1504 } else {
1505 if(origownvtx) {
1506 delete origownvtx; origownvtx=NULL;
1507 }
1508 }
1509 return;
1510}
1511//--------------------------------------------------------------------------
1512Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1513{
1514 //
1515 // Checks if this candidate is matched to MC signal
1516 //
1517
1518 if(!aod) return kFALSE;
1519
1520 // get the MC array
1521 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1522
1523 if(!mcArray) return kFALSE;
1524
1525 // try to match
1526 Int_t label = d->MatchToMC(pdg,mcArray);
1527
1528 if(label>=0) {
1529 //printf("MATCH!\n");
1530 return kTRUE;
1531 }
1532
1533 return kFALSE;
1534}
1535
1536
4bdc0315 1537//--------------------------------------------------------------------------
1538Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1539 // recompute event primary vertex from AOD tracks
1540
1541 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1542 vertexer->SetITSMode();
1543 vertexer->SetMinClusters(3);
1544
1545 AliAODVertex* pvtx=event->GetPrimaryVertex();
1546 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1547 Float_t diamondcovxy[3];
1548 event->GetDiamondCovXY(diamondcovxy);
1549 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1550 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1551 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1552 vertexer->SetVtxStart(diamond);
1553 delete diamond; diamond=NULL;
1554 }
1555
1556 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1557 if(!vertexESD) return kFALSE;
1558 if(vertexESD->GetNContributors()<=0) {
1559 //AliDebug(2,"vertexing failed");
1560 delete vertexESD; vertexESD=NULL;
1561 return kFALSE;
1562 }
1563 delete vertexer; vertexer=NULL;
1564
1565 // convert to AliAODVertex
1566 Double_t pos[3],cov[6],chi2perNDF;
1567 vertexESD->GetXYZ(pos); // position
1568 vertexESD->GetCovMatrix(cov); //covariance matrix
1569 chi2perNDF = vertexESD->GetChi2toNDF();
1570 delete vertexESD; vertexESD=NULL;
1571
1572 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1573 pvtx->SetChi2perNDF(chi2perNDF);
1574 pvtx->SetCovMatrix(cov);
1575
1576 return kTRUE;
1577}