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