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