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