1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // Base class for cuts on AOD reconstructed heavy-flavour decay
20 // Author: A.Dainese, andrea.dainese@pd.infn.it
21 /////////////////////////////////////////////////////////////
22 #include <Riostream.h>
24 #include "AliVEvent.h"
25 #include "AliESDEvent.h"
26 #include "AliAODEvent.h"
27 #include "AliVVertex.h"
28 #include "AliESDVertex.h"
30 #include "AliAODVertex.h"
31 #include "AliESDtrack.h"
32 #include "AliAODTrack.h"
33 #include "AliESDtrackCuts.h"
34 #include "AliCentrality.h"
35 #include "AliAODRecoDecayHF.h"
36 #include "AliAnalysisVertexingHF.h"
37 #include "AliRDHFCuts.h"
41 //--------------------------------------------------------------------------
42 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
43 AliAnalysisCuts(name,title),
48 fMinSPDMultiplicity(0),
64 fRemoveDaughtersFromPrimary(kFALSE),
74 // Default Constructor
77 //--------------------------------------------------------------------------
78 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
79 AliAnalysisCuts(source),
80 fMinVtxType(source.fMinVtxType),
81 fMinVtxContr(source.fMinVtxContr),
82 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
83 fMaxVtxZ(source.fMaxVtxZ),
84 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
85 fTriggerMask(source.fTriggerMask),
87 fnPtBins(source.fnPtBins),
88 fnPtBinLimits(source.fnPtBinLimits),
90 fnVars(source.fnVars),
92 fnVarsForOpt(source.fnVarsForOpt),
94 fGlobalIndex(source.fGlobalIndex),
97 fUsePID(source.fUsePID),
99 fWhyRejection(source.fWhyRejection),
100 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
101 fOptPileup(source.fOptPileup),
102 fMinContrPileup(source.fMinContrPileup),
103 fMinDzPileup(source.fMinDzPileup),
104 fUseCentrality(source.fUseCentrality),
105 fMinCentrality(source.fMinCentrality),
106 fMaxCentrality(source.fMaxCentrality),
107 fFixRefs(source.fFixRefs)
112 cout<<"Copy constructor"<<endl;
113 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
114 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
115 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
116 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
117 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
118 if(source.fPidHF) SetPidHF(source.fPidHF);
122 //--------------------------------------------------------------------------
123 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
126 // assignment operator
128 if(&source == this) return *this;
130 AliAnalysisCuts::operator=(source);
132 fMinVtxType=source.fMinVtxType;
133 fMinVtxContr=source.fMinVtxContr;
134 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
135 fMaxVtxZ=source.fMaxVtxZ;
136 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
137 fTriggerMask=source.fTriggerMask;
138 fnPtBins=source.fnPtBins;
139 fnVars=source.fnVars;
140 fGlobalIndex=source.fGlobalIndex;
141 fnVarsForOpt=source.fnVarsForOpt;
142 fUsePID=source.fUsePID;
143 SetPidHF(source.GetPidHF());
144 fWhyRejection=source.fWhyRejection;
145 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
146 fOptPileup=source.fOptPileup;
147 fMinContrPileup=source.fMinContrPileup;
148 fMinDzPileup=source.fMinDzPileup;
149 fUseCentrality=source.fUseCentrality;
150 fMinCentrality=source.fMinCentrality;
151 fMaxCentrality=source.fMaxCentrality;
152 fFixRefs=source.fFixRefs;
154 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
155 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
156 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
157 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
158 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
163 //--------------------------------------------------------------------------
164 AliRDHFCuts::~AliRDHFCuts() {
166 // Default Destructor
168 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
169 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
170 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
175 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
181 //---------------------------------------------------------------------------
182 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
186 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
188 // TEMPORARY FIX FOR REFERENCES
189 // Fix references to daughter tracks
191 AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
192 fixer->FixReferences((AliAODEvent*)event);
200 // multiplicity cuts no implemented yet
202 const AliVVertex *vertex = event->GetPrimaryVertex();
204 if(!vertex) return kFALSE;
206 TString title=vertex->GetTitle();
207 if(title.Contains("Z") && fMinVtxType>1) return kFALSE;
208 if(title.Contains("3D") && fMinVtxType>2) return kFALSE;
210 if(vertex->GetNContributors()<fMinVtxContr) return kFALSE;
212 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) return kFALSE;
214 // switch to settings for 1-pad cls in TPC
216 if(event->GetRunNumber()>121693 && event->GetRunNumber()<136851)
217 fPidHF->SetOnePad(kTRUE);
218 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517)
219 fPidHF->SetPbPb(kTRUE);
222 if(fOptPileup==kRejectPileupEvent){
223 Int_t cutc=(Int_t)fMinContrPileup;
224 Double_t cutz=(Double_t)fMinDzPileup;
225 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
231 //centrality selection
232 if (!(fUseCentrality==kCentOff)){
233 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
234 AliWarning("Centrality estimator not valid");
238 Float_t centvalue=GetCentrality((AliAODEvent*)event); if (centvalue<0.){
239 if (fWhyRejection==3) return kFALSE;
244 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
255 //---------------------------------------------------------------------------
256 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
258 // Daughter tracks selection
260 if(!fTrackCuts) return kTRUE;
262 Int_t ndaughters = d->GetNDaughters();
263 AliAODVertex *vAOD = d->GetPrimaryVtx();
264 Double_t pos[3],cov[6];
266 vAOD->GetCovarianceMatrix(cov);
267 const AliESDVertex vESD(pos,cov,100.,100);
271 for(Int_t idg=0; idg<ndaughters; idg++) {
272 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
273 if(!dgTrack) {retval = kFALSE; continue;}
274 //printf("charge %d\n",dgTrack->Charge());
275 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
277 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
282 //---------------------------------------------------------------------------
283 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
285 // Convert to ESDtrack, relate to vertex and check cuts
287 if(!cuts) return kTRUE;
291 // convert to ESD track here
292 AliESDtrack esdTrack(track);
293 // needed to calculate the impact parameters
294 esdTrack.RelateToVertex(primary,0.,3.);
295 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
297 if(fOptPileup==kRejectTracksFromPileupVertex){
299 // we need either to have here the AOD Event,
300 // or to have the pileup vertex object
304 //---------------------------------------------------------------------------
305 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
309 delete [] fPtBinLimits;
311 printf("Changing the pt bins\n");
314 if(nPtBinLimits != fnPtBins+1){
315 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
316 SetNPtBins(nPtBinLimits-1);
319 fnPtBinLimits = nPtBinLimits;
321 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
322 fPtBinLimits = new Float_t[fnPtBinLimits];
323 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
327 //---------------------------------------------------------------------------
328 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
329 // Set the variable names
334 //printf("Changing the variable names\n");
337 printf("Wrong number of variables: it has to be %d\n",fnVars);
341 fVarNames = new TString[nVars];
342 fIsUpperCut = new Bool_t[nVars];
343 for(Int_t iv=0; iv<nVars; iv++) {
344 fVarNames[iv] = varNames[iv];
345 fIsUpperCut[iv] = isUpperCut[iv];
350 //---------------------------------------------------------------------------
351 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
352 // Set the variables to be used for cuts optimization
355 delete [] fVarsForOpt;
357 //printf("Changing the variables for cut optimization\n");
360 if(nVars==0){//!=fnVars) {
361 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
366 fVarsForOpt = new Bool_t[fnVars];
367 for(Int_t iv=0; iv<fnVars; iv++) {
368 fVarsForOpt[iv]=forOpt[iv];
369 if(fVarsForOpt[iv]) fnVarsForOpt++;
375 //---------------------------------------------------------------------------
376 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
378 // set centrality estimator
381 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
387 //---------------------------------------------------------------------------
388 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
393 printf("Wrong number of variables: it has to be %d\n",fnVars);
396 if(nPtBins!=fnPtBins) {
397 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
401 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
404 for(Int_t iv=0; iv<fnVars; iv++) {
406 for(Int_t ib=0; ib<fnPtBins; ib++) {
409 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
410 cout<<"Overflow, exit..."<<endl;
414 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
420 //---------------------------------------------------------------------------
421 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
425 if(glIndex != fGlobalIndex){
426 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
429 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
431 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
432 fCutsRD[iGl] = cutsRDGlob[iGl];
436 //---------------------------------------------------------------------------
437 void AliRDHFCuts::PrintAll() const {
439 // print all cuts values
442 printf("Minimum vtx type %d\n",fMinVtxType);
443 printf("Minimum vtx contr %d\n",fMinVtxContr);
444 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
445 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
446 printf("Use PID %d\n",(Int_t)fUsePID);
447 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
448 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
449 if(fOptPileup==1) printf(" -- Reject pileup event");
450 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
451 if(fUseCentrality>0) {
452 TString estimator="";
453 if(fUseCentrality==1) estimator = "V0";
454 if(fUseCentrality==2) estimator = "Tracks";
455 if(fUseCentrality==3) estimator = "Tracklets";
456 if(fUseCentrality==4) estimator = "SPD clusters outer";
457 printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
461 cout<<"Array of variables"<<endl;
462 for(Int_t iv=0;iv<fnVars;iv++){
463 cout<<fVarNames[iv]<<"\t";
468 cout<<"Array of optimization"<<endl;
469 for(Int_t iv=0;iv<fnVars;iv++){
470 cout<<fVarsForOpt[iv]<<"\t";
475 cout<<"Array of upper/lower cut"<<endl;
476 for(Int_t iv=0;iv<fnVars;iv++){
477 cout<<fIsUpperCut[iv]<<"\t";
482 cout<<"Array of ptbin limits"<<endl;
483 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
484 cout<<fPtBinLimits[ib]<<"\t";
489 cout<<"Matrix of cuts"<<endl;
490 for(Int_t iv=0;iv<fnVars;iv++){
491 for(Int_t ib=0;ib<fnPtBins;ib++){
492 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
500 //---------------------------------------------------------------------------
501 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
506 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
511 //cout<<"Initialization..."<<endl;
512 cutsRD=new Float_t*[fnVars];
513 for(iv=0; iv<fnVars; iv++) {
514 cutsRD[iv] = new Float_t[fnPtBins];
518 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
519 GetVarPtIndex(iGlobal,iv,ib);
520 cutsRD[iv][ib] = fCutsRD[iGlobal];
526 //---------------------------------------------------------------------------
527 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
529 // give the global index from variable and pt bin
531 return iPtBin*fnVars+iVar;
534 //---------------------------------------------------------------------------
535 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
537 //give the index of the variable and of the pt bin from the global index
539 iPtBin=(Int_t)iGlob/fnVars;
545 //---------------------------------------------------------------------------
546 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
548 //give the pt bin where the pt lies.
551 for (Int_t i=0;i<fnPtBins;i++){
552 if(pt<fPtBinLimits[i+1]) {
559 //-------------------------------------------------------------------
560 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
562 // Give the value of cut set for the variable iVar and the pt bin iPtBin
565 cout<<"Cuts not iniziaisez yet"<<endl;
568 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
570 //-------------------------------------------------------------------
571 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) const {
573 // Get centrality percentile
575 AliAODHeader *header=aodEvent->GetHeader();
576 AliCentrality *centrality=header->GetCentralityP();
578 if(!centrality) return cent;
580 if (estimator==kCentV0M) cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
582 if (estimator==kCentTRK) cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
584 if (estimator==kCentTKL) cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
586 if (estimator==kCentCL1) cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
588 AliWarning("Centrality estimator not valid");
597 //-------------------------------------------------------------------
598 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
600 // Compare two cuts objects
603 Bool_t areEqual=kTRUE;
605 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
607 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
609 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
611 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
613 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
615 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
617 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
619 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
621 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
623 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;}
627 for(Int_t iv=0;iv<fnVars;iv++) {
628 for(Int_t ib=0;ib<fnPtBins;ib++) {
629 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
630 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
639 //---------------------------------------------------------------------------
640 void AliRDHFCuts::MakeTable() const {
642 // print cuts values in table format
645 TString ptString = "pT range";
646 if(fVarNames && fPtBinLimits && fCutsRD){
647 TString firstLine(Form("* %-15s",ptString.Data()));
648 for (Int_t ivar=0; ivar<fnVars; ivar++){
649 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
654 Printf("%s",firstLine.Data());
656 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
658 if (ipt==fnPtBins-1){
659 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
662 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
664 for (Int_t ivar=0; ivar<fnVars; ivar++){
665 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
667 Printf("%s",line.Data());