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