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"
29 #include "AliAODVertex.h"
30 #include "AliESDtrack.h"
31 #include "AliAODTrack.h"
32 #include "AliESDtrackCuts.h"
33 #include "AliAODRecoDecayHF.h"
34 #include "AliRDHFCuts.h"
38 //--------------------------------------------------------------------------
39 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
40 AliAnalysisCuts(name,title),
44 fMinSPDMultiplicity(0),
60 fRemoveDaughtersFromPrimary(kFALSE),
66 // Default Constructor
69 //--------------------------------------------------------------------------
70 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
71 AliAnalysisCuts(source),
72 fMinVtxType(source.fMinVtxType),
73 fMinVtxContr(source.fMinVtxContr),
74 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
75 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
76 fTriggerMask(source.fTriggerMask),
78 fnPtBins(source.fnPtBins),
79 fnPtBinLimits(source.fnPtBinLimits),
81 fnVars(source.fnVars),
83 fnVarsForOpt(source.fnVarsForOpt),
85 fGlobalIndex(source.fGlobalIndex),
88 fUsePID(source.fUsePID),
90 fWhyRejection(source.fWhyRejection),
91 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
92 fOptPileup(source.fOptPileup),
93 fMinContrPileup(source.fMinContrPileup),
94 fMinDzPileup(source.fMinDzPileup)
99 cout<<"Copy constructor"<<endl;
100 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
101 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
102 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
103 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
104 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
105 if(source.fPidHF) SetPidHF(source.fPidHF);
109 //--------------------------------------------------------------------------
110 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
113 // assignment operator
115 if(&source == this) return *this;
117 AliAnalysisCuts::operator=(source);
119 fMinVtxType=source.fMinVtxType;
120 fMinVtxContr=source.fMinVtxContr;
121 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
122 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
123 fTriggerMask=source.fTriggerMask;
124 fnPtBins=source.fnPtBins;
125 fnVars=source.fnVars;
126 fGlobalIndex=source.fGlobalIndex;
127 fnVarsForOpt=source.fnVarsForOpt;
128 fUsePID=source.fUsePID;
129 SetPidHF(source.GetPidHF());
130 fWhyRejection=source.fWhyRejection;
131 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
132 fOptPileup=source.fOptPileup;
133 fMinContrPileup=source.fMinContrPileup;
134 fMinDzPileup=source.fMinDzPileup;
136 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
137 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
138 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
139 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
140 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
145 //--------------------------------------------------------------------------
146 AliRDHFCuts::~AliRDHFCuts() {
148 // Default Destructor
150 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
151 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
152 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
157 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
163 //---------------------------------------------------------------------------
164 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
168 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
172 // multiplicity cuts no implemented yet
176 const AliVVertex *vertex = event->GetPrimaryVertex();
178 if(!vertex) return kFALSE;
180 TString title=vertex->GetTitle();
181 if(title.Contains("Z") && fMinVtxType>1) return kFALSE;
182 if(title.Contains("3D") && fMinVtxType>2) return kFALSE;
184 if(vertex->GetNContributors()<fMinVtxContr) return kFALSE;
185 if(fOptPileup==kRejectPileupEvent){
186 Int_t cutc=(Int_t)fMinContrPileup;
187 Double_t cutz=(Double_t)fMinDzPileup;
188 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
196 //---------------------------------------------------------------------------
197 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
199 // Daughter tracks selection
201 if(!fTrackCuts) return kTRUE;
203 Int_t ndaughters = d->GetNDaughters();
204 AliAODVertex *vAOD = d->GetPrimaryVtx();
205 Double_t pos[3],cov[6];
207 vAOD->GetCovarianceMatrix(cov);
208 const AliESDVertex vESD(pos,cov,100.,100);
212 for(Int_t idg=0; idg<ndaughters; idg++) {
213 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
214 if(!dgTrack) retval = kFALSE;
215 //printf("charge %d\n",dgTrack->Charge());
216 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
218 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
223 //---------------------------------------------------------------------------
224 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
226 // Convert to ESDtrack, relate to vertex and check cuts
228 if(!cuts) return kTRUE;
232 // convert to ESD track here
233 AliESDtrack esdTrack(track);
234 // needed to calculate the impact parameters
235 esdTrack.RelateToVertex(primary,0.,3.);
236 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
238 if(fOptPileup==kRejectTracksFromPileupVertex){
240 // we need either to have here the AOD Event,
241 // or to have the pileup vertex object
245 //---------------------------------------------------------------------------
246 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
250 delete [] fPtBinLimits;
252 printf("Changing the pt bins\n");
255 if(nPtBinLimits != fnPtBins+1){
256 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
257 SetNPtBins(nPtBinLimits-1);
260 fnPtBinLimits = nPtBinLimits;
262 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
263 fPtBinLimits = new Float_t[fnPtBinLimits];
264 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
268 //---------------------------------------------------------------------------
269 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
270 // Set the variable names
275 //printf("Changing the variable names\n");
278 printf("Wrong number of variables: it has to be %d\n",fnVars);
282 fVarNames = new TString[nVars];
283 fIsUpperCut = new Bool_t[nVars];
284 for(Int_t iv=0; iv<nVars; iv++) {
285 fVarNames[iv] = varNames[iv];
286 fIsUpperCut[iv] = isUpperCut[iv];
291 //---------------------------------------------------------------------------
292 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
293 // Set the variables to be used for cuts optimization
296 delete [] fVarsForOpt;
298 //printf("Changing the variables for cut optimization\n");
301 if(nVars==0){//!=fnVars) {
302 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
307 fVarsForOpt = new Bool_t[fnVars];
308 for(Int_t iv=0; iv<fnVars; iv++) {
309 fVarsForOpt[iv]=forOpt[iv];
310 if(fVarsForOpt[iv]) fnVarsForOpt++;
315 //---------------------------------------------------------------------------
316 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
321 printf("Wrong number of variables: it has to be %d\n",fnVars);
324 if(nPtBins!=fnPtBins) {
325 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
329 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
332 for(Int_t iv=0; iv<fnVars; iv++) {
334 for(Int_t ib=0; ib<fnPtBins; ib++) {
337 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
338 cout<<"Overflow, exit..."<<endl;
342 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
348 //---------------------------------------------------------------------------
349 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
353 if(glIndex != fGlobalIndex){
354 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
357 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
359 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
360 fCutsRD[iGl] = cutsRDGlob[iGl];
364 //---------------------------------------------------------------------------
365 void AliRDHFCuts::PrintAll() const {
367 // print all cuts values
370 printf("Minimum vtx type %d\n",fMinVtxType);
371 printf("Minimum vtx contr %d\n",fMinVtxContr);
372 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
373 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
374 printf("Use PID %d\n",(Int_t)fUsePID);
375 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
378 cout<<"Array of variables"<<endl;
379 for(Int_t iv=0;iv<fnVars;iv++){
380 cout<<fVarNames[iv]<<"\t";
385 cout<<"Array of optimization"<<endl;
386 for(Int_t iv=0;iv<fnVars;iv++){
387 cout<<fVarsForOpt[iv]<<"\t";
392 cout<<"Array of upper/lower cut"<<endl;
393 for(Int_t iv=0;iv<fnVars;iv++){
394 cout<<fIsUpperCut[iv]<<"\t";
399 cout<<"Array of ptbin limits"<<endl;
400 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
401 cout<<fPtBinLimits[ib]<<"\t";
406 cout<<"Matrix of cuts"<<endl;
407 for(Int_t iv=0;iv<fnVars;iv++){
408 for(Int_t ib=0;ib<fnPtBins;ib++){
409 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
417 //---------------------------------------------------------------------------
418 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
423 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
428 //cout<<"Initialization..."<<endl;
429 cutsRD=new Float_t*[fnVars];
430 for(iv=0; iv<fnVars; iv++) {
431 cutsRD[iv] = new Float_t[fnPtBins];
435 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
436 GetVarPtIndex(iGlobal,iv,ib);
437 cutsRD[iv][ib] = fCutsRD[iGlobal];
443 //---------------------------------------------------------------------------
444 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
446 // give the global index from variable and pt bin
448 return iPtBin*fnVars+iVar;
451 //---------------------------------------------------------------------------
452 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
454 //give the index of the variable and of the pt bin from the global index
456 iPtBin=(Int_t)iGlob/fnVars;
462 //---------------------------------------------------------------------------
463 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
465 //give the pt bin where the pt lies.
468 for (Int_t i=0;i<fnPtBins;i++){
469 if(pt<fPtBinLimits[i+1]) {
476 //-------------------------------------------------------------------
477 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
479 // Give the value of cut set for the variable iVar and the pt bin iPtBin
482 cout<<"Cuts not iniziaisez yet"<<endl;
485 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
487 //-------------------------------------------------------------------
488 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
490 // Compare two cuts objects
493 Bool_t areEqual=kTRUE;
495 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
497 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
499 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
501 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
503 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
505 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
507 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
509 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
511 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
513 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;}
517 for(Int_t iv=0;iv<fnVars;iv++) {
518 for(Int_t ib=0;ib<fnPtBins;ib++) {
519 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
520 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";