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