]> 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
Possibility to recompute primary vertex
[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"
46
47ClassImp(AliRDHFCuts)
48
49//--------------------------------------------------------------------------
50AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
51AliAnalysisCuts(name,title),
52fMinVtxType(3),
53fMinVtxContr(1),
54fMaxVtxRedChi2(1e6),
55fMaxVtxZ(10.),
56fMinSPDMultiplicity(0),
57fTriggerMask(AliVEvent::kAnyINT),
58fUseOnlyOneTrigger(kFALSE),
59fTriggerClass("CINT1"),
60fTrackCuts(0),
61fnPtBins(1),
62fnPtBinLimits(1),
63fPtBinLimits(0),
64fnVars(1),
65fVarNames(0),
66fnVarsForOpt(0),
67fVarsForOpt(0),
68fGlobalIndex(1),
69fCutsRD(0),
70fIsUpperCut(0),
71fUsePID(kFALSE),
72fUseAOD049(kFALSE),
73fPidHF(0),
74fWhyRejection(0),
75fEvRejectionBits(0),
76fRemoveDaughtersFromPrimary(kFALSE),
4bdc0315 77fRecomputePrimVertex(kFALSE),
0c22e2ac 78fUseMCVertex(kFALSE),
79fUsePhysicsSelection(kTRUE),
80fOptPileup(0),
81fMinContrPileup(3),
82fMinDzPileup(0.6),
83fUseCentrality(0),
84fMinCentrality(0.),
85fMaxCentrality(100.),
86fFixRefs(kFALSE),
87fIsSelectedCuts(0),
88fIsSelectedPID(0),
89fMinPtCand(-1.),
90fMaxPtCand(100000.),
91fKeepSignalMC(kFALSE),
92fIsCandTrackSPDFirst(kFALSE),
93fMaxPtCandTrackSPDFirst(0.)
94{
95 //
96 // Default Constructor
97 //
98}
99//--------------------------------------------------------------------------
100AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
101 AliAnalysisCuts(source),
102 fMinVtxType(source.fMinVtxType),
103 fMinVtxContr(source.fMinVtxContr),
104 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
105 fMaxVtxZ(source.fMaxVtxZ),
106 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
107 fTriggerMask(source.fTriggerMask),
108 fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
109 fTriggerClass(source.fTriggerClass),
110 fTrackCuts(0),
111 fnPtBins(source.fnPtBins),
112 fnPtBinLimits(source.fnPtBinLimits),
113 fPtBinLimits(0),
114 fnVars(source.fnVars),
115 fVarNames(0),
116 fnVarsForOpt(source.fnVarsForOpt),
117 fVarsForOpt(0),
118 fGlobalIndex(source.fGlobalIndex),
119 fCutsRD(0),
120 fIsUpperCut(0),
121 fUsePID(source.fUsePID),
122 fUseAOD049(source.fUseAOD049),
123 fPidHF(0),
124 fWhyRejection(source.fWhyRejection),
125 fEvRejectionBits(source.fEvRejectionBits),
126 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
4bdc0315 127 fRecomputePrimVertex(source.fRecomputePrimVertex),
0c22e2ac 128 fUseMCVertex(source.fUseMCVertex),
129 fUsePhysicsSelection(source.fUsePhysicsSelection),
130 fOptPileup(source.fOptPileup),
131 fMinContrPileup(source.fMinContrPileup),
132 fMinDzPileup(source.fMinDzPileup),
133 fUseCentrality(source.fUseCentrality),
134 fMinCentrality(source.fMinCentrality),
135 fMaxCentrality(source.fMaxCentrality),
136 fFixRefs(source.fFixRefs),
137 fIsSelectedCuts(source.fIsSelectedCuts),
138 fIsSelectedPID(source.fIsSelectedPID),
139 fMinPtCand(source.fMinPtCand),
140 fMaxPtCand(source.fMaxPtCand),
141 fKeepSignalMC(source.fKeepSignalMC),
142 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
143 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst)
144{
145 //
146 // Copy constructor
147 //
148 cout<<"Copy constructor"<<endl;
149 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
150 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
151 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
152 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
153 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
154 if(source.fPidHF) SetPidHF(source.fPidHF);
155 PrintAll();
156
157}
158//--------------------------------------------------------------------------
159AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
160{
161 //
162 // assignment operator
163 //
164 if(&source == this) return *this;
165
166 AliAnalysisCuts::operator=(source);
167
168 fMinVtxType=source.fMinVtxType;
169 fMinVtxContr=source.fMinVtxContr;
170 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
171 fMaxVtxZ=source.fMaxVtxZ;
172 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
173 fTriggerMask=source.fTriggerMask;
174 fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
175 fTriggerClass=source.fTriggerClass;
176 fnPtBins=source.fnPtBins;
177 fnPtBinLimits=source.fnPtBinLimits;
178 fnVars=source.fnVars;
179 fGlobalIndex=source.fGlobalIndex;
180 fnVarsForOpt=source.fnVarsForOpt;
181 fUsePID=source.fUsePID;
182 fUseAOD049=source.fUseAOD049;
183 if(fPidHF) delete fPidHF;
184 fPidHF=new AliAODPidHF(*(source.GetPidHF()));
185 fWhyRejection=source.fWhyRejection;
186 fEvRejectionBits=source.fEvRejectionBits;
187 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
4bdc0315 188 fRecomputePrimVertex=source.fRecomputePrimVertex;
0c22e2ac 189 fUseMCVertex=source.fUseMCVertex;
190 fUsePhysicsSelection=source.fUsePhysicsSelection;
191 fOptPileup=source.fOptPileup;
192 fMinContrPileup=source.fMinContrPileup;
193 fMinDzPileup=source.fMinDzPileup;
194 fUseCentrality=source.fUseCentrality;
195 fMinCentrality=source.fMinCentrality;
196 fMaxCentrality=source.fMaxCentrality;
197 fFixRefs=source.fFixRefs;
198 fIsSelectedCuts=source.fIsSelectedCuts;
199 fIsSelectedPID=source.fIsSelectedPID;
200 fMinPtCand=source.fMinPtCand;
201 fMaxPtCand=source.fMaxPtCand;
202 fKeepSignalMC=source.fKeepSignalMC;
203 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
204 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
205
206 if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
207 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
208 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
209 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
210 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
211 PrintAll();
212
213 return *this;
214}
215//--------------------------------------------------------------------------
216AliRDHFCuts::~AliRDHFCuts() {
217 //
218 // Default Destructor
219 //
220 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
221 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
222 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
223 if(fCutsRD) {
224 delete [] fCutsRD;
225 fCutsRD=0;
226 }
227 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
228 if(fPidHF){
229 delete fPidHF;
230 fPidHF=0;
231 }
232}
233//---------------------------------------------------------------------------
234Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
235 //
236 // Centrality selection
237 //
238 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
239 AliWarning("Centrality estimator not valid");
240 return 3;
241 }else{
242 Float_t centvalue=GetCentrality((AliAODEvent*)event);
243 if (centvalue<-998.){//-999 if no centralityP
244 return 0;
245 }else{
246 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
247 return 2;
248 }
249 }
250 }
251 return 0;
252}
253//---------------------------------------------------------------------------
254Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
255 //
256 // Event selection
257 //
258 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
259
260 fWhyRejection=0;
261 fEvRejectionBits=0;
262 Bool_t accept=kTRUE;
263
4bdc0315 264 if(fRecomputePrimVertex){
265 Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
266 if(!vertOK){
267 fWhyRejection=6;
268 return kFALSE;
269 }
270 }
271
0c22e2ac 272 // check if it's MC
273 Bool_t isMC=kFALSE;
274 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
275 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
276
277 // settings for the TPC dE/dx BB parameterization
278 if(fPidHF && fPidHF->GetOldPid()) {
279 // pp, from LHC10d onwards
280 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
281 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
282 // pp, 2011 low energy run
283 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
284 fPidHF->SetppLowEn2011(kTRUE);
285 fPidHF->SetOnePad(kFALSE);
286 }
287 // PbPb LHC10h
288 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
289 // MC
290 if(isMC) fPidHF->SetMC(kTRUE);
291 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
292 fPidHF->SetMClowenpp2011(kTRUE);
293 }
294 else if(fPidHF && !fPidHF->GetOldPid()) {
295 if(fPidHF->GetPidResponse()==0x0){
296 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
297 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
298 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
299 fPidHF->SetPidResponse(pidResp);
300 }
301 }
302
303
304 // trigger class
305 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
306 // don't do for MC and for PbPb 2010 data
307 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
308 if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {
309 fWhyRejection=5;
310 fEvRejectionBits+=1<<kNotSelTrigger;
311 accept=kFALSE;
312 }
313 }
314
315 // TEMPORARY FIX FOR REFERENCES
316 // Fix references to daughter tracks
317 // if(fFixRefs) {
318 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
319 // fixer->FixReferences((AliAODEvent*)event);
320 // delete fixer;
321 // }
322 //
323
324
325 // physics selection requirements
326 if(fUsePhysicsSelection){
327 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
328 if(!isSelected) {
329 if(accept) fWhyRejection=7;
330 fEvRejectionBits+=1<<kPhysicsSelection;
331 accept=kFALSE;
332 }else{
333 if(fUseOnlyOneTrigger){
334 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
335 if(accept) fWhyRejection=7;
336 fEvRejectionBits+=1<<kPhysicsSelection;
337 accept=kFALSE;
338 }
339 }
340 }
341 }
342
343 // vertex requirements
344
345 const AliVVertex *vertex = event->GetPrimaryVertex();
346
347 if(!vertex){
348 accept=kFALSE;
349 fEvRejectionBits+=1<<kNoVertex;
350 }else{
351 TString title=vertex->GetTitle();
352 if(title.Contains("Z") && fMinVtxType>1){
353 accept=kFALSE;
354 fEvRejectionBits+=1<<kNoVertex;
355 }
356 else if(title.Contains("3D") && fMinVtxType>2){
357 accept=kFALSE;
358 fEvRejectionBits+=1<<kNoVertex;
359 }
360 if(vertex->GetNContributors()<fMinVtxContr){
361 accept=kFALSE;
362 fEvRejectionBits+=1<<kTooFewVtxContrib;
363 }
364 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
365 fEvRejectionBits+=1<<kZVtxOutFid;
366 if(accept) fWhyRejection=6;
367 accept=kFALSE;
368 }
369 }
370
371
372 // pile-up rejection
373 if(fOptPileup==kRejectPileupEvent){
374 Int_t cutc=(Int_t)fMinContrPileup;
375 Double_t cutz=(Double_t)fMinDzPileup;
376 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
377 if(accept) fWhyRejection=1;
378 fEvRejectionBits+=1<<kPileupSPD;
379 accept=kFALSE;
380 }
381 }
382
383 // centrality selection
384 if (fUseCentrality!=kCentOff) {
385 Int_t rejection=IsEventSelectedInCentrality(event);
386 if(rejection>1){
387 if(accept) fWhyRejection=rejection;
388 fEvRejectionBits+=1<<kOutsideCentrality;
389 accept=kFALSE;
390 }
391 }
392
393 return accept;
394}
395//---------------------------------------------------------------------------
396Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
397 //
398 // Daughter tracks selection
399 //
400 if(!fTrackCuts) return kTRUE;
401
402 Int_t ndaughters = d->GetNDaughters();
403 AliAODVertex *vAOD = d->GetPrimaryVtx();
404 Double_t pos[3],cov[6];
405 vAOD->GetXYZ(pos);
406 vAOD->GetCovarianceMatrix(cov);
407 const AliESDVertex vESD(pos,cov,100.,100);
408
409 Bool_t retval=kTRUE;
410
411 for(Int_t idg=0; idg<ndaughters; idg++) {
412 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
413 if(!dgTrack) {retval = kFALSE; continue;}
414 //printf("charge %d\n",dgTrack->Charge());
415 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
416
417 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
418 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
419
420 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
421 }
422
423 return retval;
424}
425//---------------------------------------------------------------------------
426Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
427 //
428 // Convert to ESDtrack, relate to vertex and check cuts
429 //
430 if(!cuts) return kTRUE;
431
432 if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
433
434 Bool_t retval=kTRUE;
435
436 // convert to ESD track here
437 AliESDtrack esdTrack(track);
438 // set the TPC cluster info
439 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
440 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
441 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
442 // needed to calculate the impact parameters
443 esdTrack.RelateToVertex(primary,0.,3.);
444
445 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
446
447 if(fOptPileup==kRejectTracksFromPileupVertex){
448 // to be implemented
449 // we need either to have here the AOD Event,
450 // or to have the pileup vertex object
451 }
452 return retval;
453}
454//---------------------------------------------------------------------------
455void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
456 // Set the pt bins
457
458 if(fPtBinLimits) {
459 delete [] fPtBinLimits;
460 fPtBinLimits = NULL;
461 printf("Changing the pt bins\n");
462 }
463
464 if(nPtBinLimits != fnPtBins+1){
465 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
466 SetNPtBins(nPtBinLimits-1);
467 }
468
469 fnPtBinLimits = nPtBinLimits;
470 SetGlobalIndex();
471 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
472 fPtBinLimits = new Float_t[fnPtBinLimits];
473 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
474
475 return;
476}
477//---------------------------------------------------------------------------
478void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
479 // Set the variable names
480
481 if(fVarNames) {
482 delete [] fVarNames;
483 fVarNames = NULL;
484 //printf("Changing the variable names\n");
485 }
486 if(nVars!=fnVars){
487 printf("Wrong number of variables: it has to be %d\n",fnVars);
488 return;
489 }
490 //fnVars=nVars;
491 fVarNames = new TString[nVars];
492 fIsUpperCut = new Bool_t[nVars];
493 for(Int_t iv=0; iv<nVars; iv++) {
494 fVarNames[iv] = varNames[iv];
495 fIsUpperCut[iv] = isUpperCut[iv];
496 }
497
498 return;
499}
500//---------------------------------------------------------------------------
501void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
502 // Set the variables to be used for cuts optimization
503
504 if(fVarsForOpt) {
505 delete [] fVarsForOpt;
506 fVarsForOpt = NULL;
507 //printf("Changing the variables for cut optimization\n");
508 }
509
510 if(nVars==0){//!=fnVars) {
511 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
512 return;
513 }
514
515 fnVarsForOpt = 0;
516 fVarsForOpt = new Bool_t[fnVars];
517 for(Int_t iv=0; iv<fnVars; iv++) {
518 fVarsForOpt[iv]=forOpt[iv];
519 if(fVarsForOpt[iv]) fnVarsForOpt++;
520 }
521
522 return;
523}
524
525//---------------------------------------------------------------------------
526void AliRDHFCuts::SetUseCentrality(Int_t flag) {
527 //
528 // set centrality estimator
529 //
530 fUseCentrality=flag;
531 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
532
533 return;
534}
535
536
537//---------------------------------------------------------------------------
538void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
539 //
540 // store the cuts
541 //
542 if(nVars!=fnVars) {
543 printf("Wrong number of variables: it has to be %d\n",fnVars);
544 AliFatal("exiting");
545 }
546 if(nPtBins!=fnPtBins) {
547 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
548 AliFatal("exiting");
549 }
550
551 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
552
553
554 for(Int_t iv=0; iv<fnVars; iv++) {
555
556 for(Int_t ib=0; ib<fnPtBins; ib++) {
557
558 //check
559 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
560 cout<<"Overflow, exit..."<<endl;
561 return;
562 }
563
564 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
565
566 }
567 }
568 return;
569}
570//---------------------------------------------------------------------------
571void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
572 //
573 // store the cuts
574 //
575 if(glIndex != fGlobalIndex){
576 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
577 AliFatal("exiting");
578 }
579 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
580
581 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
582 fCutsRD[iGl] = cutsRDGlob[iGl];
583 }
584 return;
585}
586//---------------------------------------------------------------------------
587void AliRDHFCuts::PrintAll() const {
588 //
589 // print all cuts values
590 //
591
592 printf("Minimum vtx type %d\n",fMinVtxType);
593 printf("Minimum vtx contr %d\n",fMinVtxContr);
594 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
595 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
596 printf("Use PID %d\n",(Int_t)fUsePID);
597 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
4bdc0315 598 printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
0c22e2ac 599 printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
600 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
601 if(fOptPileup==1) printf(" -- Reject pileup event");
602 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
603 if(fUseCentrality>0) {
604 TString estimator="";
605 if(fUseCentrality==1) estimator = "V0";
606 if(fUseCentrality==2) estimator = "Tracks";
607 if(fUseCentrality==3) estimator = "Tracklets";
608 if(fUseCentrality==4) estimator = "SPD clusters outer";
609 printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
610 }
611 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
612
613 if(fVarNames){
614 cout<<"Array of variables"<<endl;
615 for(Int_t iv=0;iv<fnVars;iv++){
616 cout<<fVarNames[iv]<<"\t";
617 }
618 cout<<endl;
619 }
620 if(fVarsForOpt){
621 cout<<"Array of optimization"<<endl;
622 for(Int_t iv=0;iv<fnVars;iv++){
623 cout<<fVarsForOpt[iv]<<"\t";
624 }
625 cout<<endl;
626 }
627 if(fIsUpperCut){
628 cout<<"Array of upper/lower cut"<<endl;
629 for(Int_t iv=0;iv<fnVars;iv++){
630 cout<<fIsUpperCut[iv]<<"\t";
631 }
632 cout<<endl;
633 }
634 if(fPtBinLimits){
635 cout<<"Array of ptbin limits"<<endl;
636 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
637 cout<<fPtBinLimits[ib]<<"\t";
638 }
639 cout<<endl;
640 }
641 if(fCutsRD){
642 cout<<"Matrix of cuts"<<endl;
643 for(Int_t iv=0;iv<fnVars;iv++){
644 for(Int_t ib=0;ib<fnPtBins;ib++){
645 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
646 }
647 cout<<endl;
648 }
649 cout<<endl;
650 }
651 return;
652}
653
654//--------------------------------------------------------------------------
655void AliRDHFCuts::PrintTrigger() const{
656
657 cout<<" Trigger selection pattern: ";
658 if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
659 if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
660 if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
661 if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
662 if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
663 if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
664 if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
665 if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
666 if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
667 if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
668 cout << endl<< endl;
669
670}
671
672//---------------------------------------------------------------------------
673void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
674 //
675 // get the cuts
676 //
677
678 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
679
680
681 Int_t iv,ib;
682 if(!cutsRD) {
683 //cout<<"Initialization..."<<endl;
684 cutsRD=new Float_t*[fnVars];
685 for(iv=0; iv<fnVars; iv++) {
686 cutsRD[iv] = new Float_t[fnPtBins];
687 }
688 }
689
690 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
691 GetVarPtIndex(iGlobal,iv,ib);
692 cutsRD[iv][ib] = fCutsRD[iGlobal];
693 }
694
695 return;
696}
697
698//---------------------------------------------------------------------------
699Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
700 //
701 // give the global index from variable and pt bin
702 //
703 return iPtBin*fnVars+iVar;
704}
705
706//---------------------------------------------------------------------------
707void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
708 //
709 //give the index of the variable and of the pt bin from the global index
710 //
711 iPtBin=(Int_t)iGlob/fnVars;
712 iVar=iGlob%fnVars;
713
714 return;
715}
716
717//---------------------------------------------------------------------------
718Int_t AliRDHFCuts::PtBin(Double_t pt) const {
719 //
720 //give the pt bin where the pt lies.
721 //
722 Int_t ptbin=-1;
723 if(pt<fPtBinLimits[0])return ptbin;
724 for (Int_t i=0;i<fnPtBins;i++){
725 if(pt<fPtBinLimits[i+1]) {
726 ptbin=i;
727 break;
728 }
729 }
730 return ptbin;
731}
732//-------------------------------------------------------------------
733Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
734 //
735 // Give the value of cut set for the variable iVar and the pt bin iPtBin
736 //
737 if(!fCutsRD){
738 cout<<"Cuts not iniziaisez yet"<<endl;
739 return 0;
740 }
741 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
742}
743//-------------------------------------------------------------------
744Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
745 //
746 // Get centrality percentile
747 //
748
749 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
750 if(mcArray) {fUseAOD049=kFALSE;}
751
752 AliAODHeader *header=aodEvent->GetHeader();
753 AliCentrality *centrality=header->GetCentralityP();
754 Float_t cent=-999.;
755 Bool_t isSelRun=kFALSE;
756 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
757 if(!centrality) return cent;
758 else{
759 if (estimator==kCentV0M){
760 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
761 if(cent<0){
762 Int_t quality = centrality->GetQuality();
763 if(quality<=1){
764 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
765 }else{
766 Int_t runnum=aodEvent->GetRunNumber();
767 for(Int_t ir=0;ir<5;ir++){
768 if(runnum==selRun[ir]){
769 isSelRun=kTRUE;
770 break;
771 }
772 }
773 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
774 }
775 }
776
777 //temporary fix for AOD049 outliers
778 if(fUseAOD049&&cent>=0){
779 Float_t v0=0;
780 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
781 v0+=aodV0->GetMTotV0A();
782 v0+=aodV0->GetMTotV0C();
783 if(cent==0&&v0<19500)return -1;//filtering issue
784 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
785 Float_t val= 1.30552 + 0.147931 * v0;
786 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};
787 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
788 }
789 }
790 else {
791 if (estimator==kCentTRK) {
792 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
793 if(cent<0){
794 Int_t quality = centrality->GetQuality();
795 if(quality<=1){
796 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
797 }else{
798 Int_t runnum=aodEvent->GetRunNumber();
799 for(Int_t ir=0;ir<5;ir++){
800 if(runnum==selRun[ir]){
801 isSelRun=kTRUE;
802 break;
803 }
804 }
805 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
806 }
807 }
808 }
809 else{
810 if (estimator==kCentTKL){
811 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
812 if(cent<0){
813 Int_t quality = centrality->GetQuality();
814 if(quality<=1){
815 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
816 }else{
817 Int_t runnum=aodEvent->GetRunNumber();
818 for(Int_t ir=0;ir<5;ir++){
819 if(runnum==selRun[ir]){
820 isSelRun=kTRUE;
821 break;
822 }
823 }
824 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
825 }
826 }
827 }
828 else{
829 if (estimator==kCentCL1){
830 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
831 if(cent<0){
832 Int_t quality = centrality->GetQuality();
833 if(quality<=1){
834 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
835 }else{
836 Int_t runnum=aodEvent->GetRunNumber();
837 for(Int_t ir=0;ir<5;ir++){
838 if(runnum==selRun[ir]){
839 isSelRun=kTRUE;
840 break;
841 }
842 }
843 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
844 }
845 }
846 }
847 else {
848 AliWarning("Centrality estimator not valid");
849
850 }
851 }
852 }
853 }
854 }
855 return cent;
856}
857//-------------------------------------------------------------------
858Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
859 //
860 // Compare two cuts objects
861 //
862
863 Bool_t areEqual=kTRUE;
864
865 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
866
867 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
868
869 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
870
871 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
872
873 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
874
875 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
876 if(fTrackCuts){
877 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
878
879 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
880
881 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
882
883 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;}
884 }
885
886 if(fCutsRD) {
887 for(Int_t iv=0;iv<fnVars;iv++) {
888 for(Int_t ib=0;ib<fnPtBins;ib++) {
889 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
890 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
891 areEqual=kFALSE;
892 }
893 }
894 }
895 }
896
897 return areEqual;
898}
899//---------------------------------------------------------------------------
900void AliRDHFCuts::MakeTable() const {
901 //
902 // print cuts values in table format
903 //
904
905 TString ptString = "pT range";
906 if(fVarNames && fPtBinLimits && fCutsRD){
907 TString firstLine(Form("* %-15s",ptString.Data()));
908 for (Int_t ivar=0; ivar<fnVars; ivar++){
909 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
910 if (ivar == fnVars){
911 firstLine+="*\n";
912 }
913 }
914 Printf("%s",firstLine.Data());
915
916 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
917 TString line;
918 if (ipt==fnPtBins-1){
919 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
920 }
921 else{
922 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
923 }
924 for (Int_t ivar=0; ivar<fnVars; ivar++){
925 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
926 }
927 Printf("%s",line.Data());
928 }
929
930 }
931
932
933 return;
934}
935//--------------------------------------------------------------------------
936Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
937 AliAODEvent *aod) const
938{
939 //
940 // Recalculate primary vertex without daughters
941 //
942
943 if(!aod) {
944 AliError("Can not remove daughters from vertex without AOD event");
945 return 0;
946 }
947
948 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
949 if(!recvtx){
950 AliDebug(2,"Removal of daughter tracks failed");
951 return kFALSE;
952 }
953
954
955 //set recalculed primary vertex
956 d->SetOwnPrimaryVtx(recvtx);
957 delete recvtx;
958
959 return kTRUE;
960}
961//--------------------------------------------------------------------------
962Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
963{
964 //
965 // Recalculate primary vertex without daughters
966 //
967
968 if(!aod) {
969 AliError("Can not get MC vertex without AOD event");
970 return kFALSE;
971 }
972
973 // load MC header
974 AliAODMCHeader *mcHeader =
975 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
976 if(!mcHeader) {
977 AliError("Can not get MC vertex without AODMCHeader event");
978 return kFALSE;
979 }
980 Double_t pos[3];
981 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
982 mcHeader->GetVertex(pos);
983 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
984
985 if(!recvtx){
986 AliDebug(2,"Removal of daughter tracks failed");
987 return kFALSE;
988 }
989
990 //set recalculed primary vertex
991 d->SetOwnPrimaryVtx(recvtx);
992
993 d->RecalculateImpPars(recvtx,aod);
994
995 delete recvtx;
996
997 return kTRUE;
998}
999//--------------------------------------------------------------------------
1000void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1001 AliAODEvent *aod,
1002 AliAODVertex *origownvtx) const
1003{
1004 //
1005 // Clean-up own primary vertex if needed
1006 //
1007
1008 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1009 d->UnsetOwnPrimaryVtx();
1010 if(origownvtx) {
1011 d->SetOwnPrimaryVtx(origownvtx);
1012 delete origownvtx; origownvtx=NULL;
1013 }
1014 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1015 } else {
1016 if(origownvtx) {
1017 delete origownvtx; origownvtx=NULL;
1018 }
1019 }
1020 return;
1021}
1022//--------------------------------------------------------------------------
1023Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1024{
1025 //
1026 // Checks if this candidate is matched to MC signal
1027 //
1028
1029 if(!aod) return kFALSE;
1030
1031 // get the MC array
1032 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1033
1034 if(!mcArray) return kFALSE;
1035
1036 // try to match
1037 Int_t label = d->MatchToMC(pdg,mcArray);
1038
1039 if(label>=0) {
1040 //printf("MATCH!\n");
1041 return kTRUE;
1042 }
1043
1044 return kFALSE;
1045}
1046
1047
4bdc0315 1048//--------------------------------------------------------------------------
1049Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1050 // recompute event primary vertex from AOD tracks
1051
1052 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1053 vertexer->SetITSMode();
1054 vertexer->SetMinClusters(3);
1055
1056 AliAODVertex* pvtx=event->GetPrimaryVertex();
1057 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1058 Float_t diamondcovxy[3];
1059 event->GetDiamondCovXY(diamondcovxy);
1060 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1061 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1062 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1063 vertexer->SetVtxStart(diamond);
1064 delete diamond; diamond=NULL;
1065 }
1066
1067 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1068 if(!vertexESD) return kFALSE;
1069 if(vertexESD->GetNContributors()<=0) {
1070 //AliDebug(2,"vertexing failed");
1071 delete vertexESD; vertexESD=NULL;
1072 return kFALSE;
1073 }
1074 delete vertexer; vertexer=NULL;
1075
1076 // convert to AliAODVertex
1077 Double_t pos[3],cov[6],chi2perNDF;
1078 vertexESD->GetXYZ(pos); // position
1079 vertexESD->GetCovMatrix(cov); //covariance matrix
1080 chi2perNDF = vertexESD->GetChi2toNDF();
1081 delete vertexESD; vertexESD=NULL;
1082
1083 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1084 pvtx->SetChi2perNDF(chi2perNDF);
1085 pvtx->SetCovMatrix(cov);
1086
1087 return kTRUE;
1088}