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