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