1 /**************************************************************************
\r
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
18 /////////////////////////////////////////////////////////////
\r
20 // Base class for cuts on AOD reconstructed heavy-flavour decay
\r
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
\r
23 /////////////////////////////////////////////////////////////
\r
24 #include <Riostream.h>
\r
26 #include "AliVEvent.h"
\r
27 #include "AliESDEvent.h"
\r
28 #include "AliAODEvent.h"
\r
29 #include "AliVVertex.h"
\r
30 #include "AliESDVertex.h"
\r
32 #include "AliAODVertex.h"
\r
33 #include "AliESDtrack.h"
\r
34 #include "AliAODTrack.h"
\r
35 #include "AliESDtrackCuts.h"
\r
36 #include "AliCentrality.h"
\r
37 #include "AliAODRecoDecayHF.h"
\r
38 #include "AliAnalysisVertexingHF.h"
\r
39 #include "AliAODMCHeader.h"
\r
40 #include "AliAODMCParticle.h"
\r
41 #include "AliRDHFCuts.h"
\r
42 #include "AliAnalysisManager.h"
\r
43 #include "AliInputEventHandler.h"
\r
44 #include "AliPIDResponse.h"
\r
46 ClassImp(AliRDHFCuts)
\r
48 //--------------------------------------------------------------------------
\r
49 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
\r
50 AliAnalysisCuts(name,title),
\r
53 fMaxVtxRedChi2(1e6),
\r
55 fMinSPDMultiplicity(0),
\r
57 fTriggerClass("CINT1"),
\r
73 fEvRejectionBits(0),
\r
74 fRemoveDaughtersFromPrimary(kFALSE),
\r
75 fUseMCVertex(kFALSE),
\r
76 fUsePhysicsSelection(kTRUE),
\r
82 fMaxCentrality(100.),
\r
87 fMaxPtCand(100000.),
\r
88 fKeepSignalMC(kFALSE),
\r
89 fIsCandTrackSPDFirst(kFALSE),
\r
90 fMaxPtCandTrackSPDFirst(0.)
\r
93 // Default Constructor
\r
96 //--------------------------------------------------------------------------
\r
97 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
\r
98 AliAnalysisCuts(source),
\r
99 fMinVtxType(source.fMinVtxType),
\r
100 fMinVtxContr(source.fMinVtxContr),
\r
101 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
\r
102 fMaxVtxZ(source.fMaxVtxZ),
\r
103 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
\r
104 fTriggerMask(source.fTriggerMask),
\r
105 fTriggerClass(source.fTriggerClass),
\r
107 fnPtBins(source.fnPtBins),
\r
108 fnPtBinLimits(source.fnPtBinLimits),
\r
110 fnVars(source.fnVars),
\r
112 fnVarsForOpt(source.fnVarsForOpt),
\r
114 fGlobalIndex(source.fGlobalIndex),
\r
117 fUsePID(source.fUsePID),
\r
118 fUseAOD049(source.fUseAOD049),
\r
120 fWhyRejection(source.fWhyRejection),
\r
121 fEvRejectionBits(source.fEvRejectionBits),
\r
122 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
\r
123 fUseMCVertex(source.fUseMCVertex),
\r
124 fUsePhysicsSelection(source.fUsePhysicsSelection),
\r
125 fOptPileup(source.fOptPileup),
\r
126 fMinContrPileup(source.fMinContrPileup),
\r
127 fMinDzPileup(source.fMinDzPileup),
\r
128 fUseCentrality(source.fUseCentrality),
\r
129 fMinCentrality(source.fMinCentrality),
\r
130 fMaxCentrality(source.fMaxCentrality),
\r
131 fFixRefs(source.fFixRefs),
\r
132 fIsSelectedCuts(source.fIsSelectedCuts),
\r
133 fIsSelectedPID(source.fIsSelectedPID),
\r
134 fMinPtCand(source.fMinPtCand),
\r
135 fMaxPtCand(source.fMaxPtCand),
\r
136 fKeepSignalMC(source.fKeepSignalMC),
\r
137 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
\r
138 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst)
\r
141 // Copy constructor
\r
143 cout<<"Copy constructor"<<endl;
\r
144 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
\r
145 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
\r
146 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
\r
147 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
\r
148 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
\r
149 if(source.fPidHF) SetPidHF(source.fPidHF);
\r
153 //--------------------------------------------------------------------------
\r
154 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
\r
157 // assignment operator
\r
159 if(&source == this) return *this;
\r
161 AliAnalysisCuts::operator=(source);
\r
163 fMinVtxType=source.fMinVtxType;
\r
164 fMinVtxContr=source.fMinVtxContr;
\r
165 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
\r
166 fMaxVtxZ=source.fMaxVtxZ;
\r
167 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
\r
168 fTriggerMask=source.fTriggerMask;
\r
169 fTriggerClass=source.fTriggerClass;
\r
170 fnPtBins=source.fnPtBins;
\r
171 fnVars=source.fnVars;
\r
172 fGlobalIndex=source.fGlobalIndex;
\r
173 fnVarsForOpt=source.fnVarsForOpt;
\r
174 fUsePID=source.fUsePID;
\r
175 fUseAOD049=source.fUseAOD049;
\r
176 SetPidHF(source.GetPidHF());
\r
177 fWhyRejection=source.fWhyRejection;
\r
178 fEvRejectionBits=source.fEvRejectionBits;
\r
179 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
\r
180 fUseMCVertex=source.fUseMCVertex;
\r
181 fUsePhysicsSelection=source.fUsePhysicsSelection;
\r
182 fOptPileup=source.fOptPileup;
\r
183 fMinContrPileup=source.fMinContrPileup;
\r
184 fMinDzPileup=source.fMinDzPileup;
\r
185 fUseCentrality=source.fUseCentrality;
\r
186 fMinCentrality=source.fMinCentrality;
\r
187 fMaxCentrality=source.fMaxCentrality;
\r
188 fFixRefs=source.fFixRefs;
\r
189 fIsSelectedCuts=source.fIsSelectedCuts;
\r
190 fIsSelectedPID=source.fIsSelectedPID;
\r
191 fMinPtCand=source.fMinPtCand;
\r
192 fMaxPtCand=source.fMaxPtCand;
\r
193 fKeepSignalMC=source.fKeepSignalMC;
\r
194 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
\r
195 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
\r
197 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
\r
198 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
\r
199 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
\r
200 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
\r
201 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
\r
206 //--------------------------------------------------------------------------
\r
207 AliRDHFCuts::~AliRDHFCuts() {
\r
209 // Default Destructor
\r
211 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
\r
212 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
\r
213 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
\r
218 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
\r
224 //---------------------------------------------------------------------------
\r
225 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
\r
227 // Centrality selection
\r
229 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
\r
230 AliWarning("Centrality estimator not valid");
\r
233 Float_t centvalue=GetCentrality((AliAODEvent*)event);
\r
234 if (centvalue<-998.){//-999 if no centralityP
\r
237 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
\r
244 //---------------------------------------------------------------------------
\r
245 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
\r
249 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
\r
252 fEvRejectionBits=0;
\r
253 Bool_t accept=kTRUE;
\r
255 // check if it's MC
\r
256 Bool_t isMC=kFALSE;
\r
257 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
\r
258 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
\r
260 // settings for the TPC dE/dx BB parameterization
\r
261 if(fPidHF && fPidHF->GetOldPid()) {
\r
262 // pp, from LHC10d onwards
\r
263 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
\r
264 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
\r
265 // pp, 2011 low energy run
\r
266 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
\r
267 fPidHF->SetppLowEn2011(kTRUE);
\r
268 fPidHF->SetOnePad(kFALSE);
\r
271 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
\r
273 if(isMC) fPidHF->SetMC(kTRUE);
\r
274 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
\r
275 fPidHF->SetMClowenpp2011(kTRUE);
\r
277 else if(fPidHF && !fPidHF->GetOldPid()) {
\r
278 if(fPidHF->GetPidResponse()==0x0){
\r
279 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
\r
280 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
\r
281 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
\r
282 fPidHF->SetPidResponse(pidResp);
\r
288 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
\r
289 // don't do for MC and for PbPb 2010 data
\r
290 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
\r
291 if(!firedTriggerClasses.Contains(fTriggerClass.Data())) {
\r
293 fEvRejectionBits+=1<<kNotSelTrigger;
\r
298 // TEMPORARY FIX FOR REFERENCES
\r
299 // Fix references to daughter tracks
\r
301 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
\r
302 // fixer->FixReferences((AliAODEvent*)event);
\r
308 // physics selection requirements
\r
309 if(fUsePhysicsSelection){
\r
310 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
\r
312 if(accept) fWhyRejection=7;
\r
313 fEvRejectionBits+=1<<kPhysicsSelection;
\r
318 // vertex requirements
\r
320 const AliVVertex *vertex = event->GetPrimaryVertex();
\r
324 fEvRejectionBits+=1<<kNoVertex;
\r
326 TString title=vertex->GetTitle();
\r
327 if(title.Contains("Z") && fMinVtxType>1){
\r
329 fEvRejectionBits+=1<<kNoVertex;
\r
331 else if(title.Contains("3D") && fMinVtxType>2){
\r
333 fEvRejectionBits+=1<<kNoVertex;
\r
335 if(vertex->GetNContributors()<fMinVtxContr){
\r
337 fEvRejectionBits+=1<<kTooFewVtxContrib;
\r
339 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
\r
340 fEvRejectionBits+=1<<kZVtxOutFid;
\r
341 if(accept) fWhyRejection=6;
\r
347 // pile-up rejection
\r
348 if(fOptPileup==kRejectPileupEvent){
\r
349 Int_t cutc=(Int_t)fMinContrPileup;
\r
350 Double_t cutz=(Double_t)fMinDzPileup;
\r
351 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
\r
352 if(accept) fWhyRejection=1;
\r
353 fEvRejectionBits+=1<<kPileupSPD;
\r
358 // centrality selection
\r
359 if (fUseCentrality!=kCentOff) {
\r
360 Int_t rejection=IsEventSelectedInCentrality(event);
\r
362 if(accept) fWhyRejection=rejection;
\r
363 fEvRejectionBits+=1<<kOutsideCentrality;
\r
370 //---------------------------------------------------------------------------
\r
371 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
\r
373 // Daughter tracks selection
\r
375 if(!fTrackCuts) return kTRUE;
\r
377 Int_t ndaughters = d->GetNDaughters();
\r
378 AliAODVertex *vAOD = d->GetPrimaryVtx();
\r
379 Double_t pos[3],cov[6];
\r
381 vAOD->GetCovarianceMatrix(cov);
\r
382 const AliESDVertex vESD(pos,cov,100.,100);
\r
384 Bool_t retval=kTRUE;
\r
386 for(Int_t idg=0; idg<ndaughters; idg++) {
\r
387 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
\r
388 if(!dgTrack) {retval = kFALSE; continue;}
\r
389 //printf("charge %d\n",dgTrack->Charge());
\r
390 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
\r
392 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
\r
393 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
\r
395 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
\r
400 //---------------------------------------------------------------------------
\r
401 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
\r
403 // Convert to ESDtrack, relate to vertex and check cuts
\r
405 if(!cuts) return kTRUE;
\r
407 Bool_t retval=kTRUE;
\r
409 // convert to ESD track here
\r
410 AliESDtrack esdTrack(track);
\r
411 // set the TPC cluster info
\r
412 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
\r
413 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
\r
414 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
\r
415 // needed to calculate the impact parameters
\r
416 esdTrack.RelateToVertex(primary,0.,3.);
\r
418 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
\r
420 if(fOptPileup==kRejectTracksFromPileupVertex){
\r
421 // to be implemented
\r
422 // we need either to have here the AOD Event,
\r
423 // or to have the pileup vertex object
\r
427 //---------------------------------------------------------------------------
\r
428 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
\r
432 delete [] fPtBinLimits;
\r
433 fPtBinLimits = NULL;
\r
434 printf("Changing the pt bins\n");
\r
437 if(nPtBinLimits != fnPtBins+1){
\r
438 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
\r
439 SetNPtBins(nPtBinLimits-1);
\r
442 fnPtBinLimits = nPtBinLimits;
\r
444 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
\r
445 fPtBinLimits = new Float_t[fnPtBinLimits];
\r
446 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
\r
450 //---------------------------------------------------------------------------
\r
451 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
\r
452 // Set the variable names
\r
455 delete [] fVarNames;
\r
457 //printf("Changing the variable names\n");
\r
460 printf("Wrong number of variables: it has to be %d\n",fnVars);
\r
464 fVarNames = new TString[nVars];
\r
465 fIsUpperCut = new Bool_t[nVars];
\r
466 for(Int_t iv=0; iv<nVars; iv++) {
\r
467 fVarNames[iv] = varNames[iv];
\r
468 fIsUpperCut[iv] = isUpperCut[iv];
\r
473 //---------------------------------------------------------------------------
\r
474 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
\r
475 // Set the variables to be used for cuts optimization
\r
478 delete [] fVarsForOpt;
\r
479 fVarsForOpt = NULL;
\r
480 //printf("Changing the variables for cut optimization\n");
\r
483 if(nVars==0){//!=fnVars) {
\r
484 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
\r
489 fVarsForOpt = new Bool_t[fnVars];
\r
490 for(Int_t iv=0; iv<fnVars; iv++) {
\r
491 fVarsForOpt[iv]=forOpt[iv];
\r
492 if(fVarsForOpt[iv]) fnVarsForOpt++;
\r
498 //---------------------------------------------------------------------------
\r
499 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
\r
501 // set centrality estimator
\r
503 fUseCentrality=flag;
\r
504 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
\r
510 //---------------------------------------------------------------------------
\r
511 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
\r
515 if(nVars!=fnVars) {
\r
516 printf("Wrong number of variables: it has to be %d\n",fnVars);
\r
517 AliFatal("exiting");
\r
519 if(nPtBins!=fnPtBins) {
\r
520 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
\r
521 AliFatal("exiting");
\r
524 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
\r
527 for(Int_t iv=0; iv<fnVars; iv++) {
\r
529 for(Int_t ib=0; ib<fnPtBins; ib++) {
\r
532 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
\r
533 cout<<"Overflow, exit..."<<endl;
\r
537 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
\r
543 //---------------------------------------------------------------------------
\r
544 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
\r
548 if(glIndex != fGlobalIndex){
\r
549 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
\r
550 AliFatal("exiting");
\r
552 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
\r
554 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
\r
555 fCutsRD[iGl] = cutsRDGlob[iGl];
\r
559 //---------------------------------------------------------------------------
\r
560 void AliRDHFCuts::PrintAll() const {
\r
562 // print all cuts values
\r
565 printf("Minimum vtx type %d\n",fMinVtxType);
\r
566 printf("Minimum vtx contr %d\n",fMinVtxContr);
\r
567 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
\r
568 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
\r
569 printf("Use PID %d\n",(Int_t)fUsePID);
\r
570 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
\r
571 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
\r
572 if(fOptPileup==1) printf(" -- Reject pileup event");
\r
573 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
\r
574 if(fUseCentrality>0) {
\r
575 TString estimator="";
\r
576 if(fUseCentrality==1) estimator = "V0";
\r
577 if(fUseCentrality==2) estimator = "Tracks";
\r
578 if(fUseCentrality==3) estimator = "Tracklets";
\r
579 if(fUseCentrality==4) estimator = "SPD clusters outer";
\r
580 printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
\r
582 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
\r
585 cout<<"Array of variables"<<endl;
\r
586 for(Int_t iv=0;iv<fnVars;iv++){
\r
587 cout<<fVarNames[iv]<<"\t";
\r
592 cout<<"Array of optimization"<<endl;
\r
593 for(Int_t iv=0;iv<fnVars;iv++){
\r
594 cout<<fVarsForOpt[iv]<<"\t";
\r
599 cout<<"Array of upper/lower cut"<<endl;
\r
600 for(Int_t iv=0;iv<fnVars;iv++){
\r
601 cout<<fIsUpperCut[iv]<<"\t";
\r
606 cout<<"Array of ptbin limits"<<endl;
\r
607 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
\r
608 cout<<fPtBinLimits[ib]<<"\t";
\r
613 cout<<"Matrix of cuts"<<endl;
\r
614 for(Int_t iv=0;iv<fnVars;iv++){
\r
615 for(Int_t ib=0;ib<fnPtBins;ib++){
\r
616 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
\r
624 //---------------------------------------------------------------------------
\r
625 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
\r
630 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
\r
635 //cout<<"Initialization..."<<endl;
\r
636 cutsRD=new Float_t*[fnVars];
\r
637 for(iv=0; iv<fnVars; iv++) {
\r
638 cutsRD[iv] = new Float_t[fnPtBins];
\r
642 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
\r
643 GetVarPtIndex(iGlobal,iv,ib);
\r
644 cutsRD[iv][ib] = fCutsRD[iGlobal];
\r
650 //---------------------------------------------------------------------------
\r
651 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
\r
653 // give the global index from variable and pt bin
\r
655 return iPtBin*fnVars+iVar;
\r
658 //---------------------------------------------------------------------------
\r
659 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
\r
661 //give the index of the variable and of the pt bin from the global index
\r
663 iPtBin=(Int_t)iGlob/fnVars;
\r
669 //---------------------------------------------------------------------------
\r
670 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
\r
672 //give the pt bin where the pt lies.
\r
675 if(pt<fPtBinLimits[0])return ptbin;
\r
676 for (Int_t i=0;i<fnPtBins;i++){
\r
677 if(pt<fPtBinLimits[i+1]) {
\r
684 //-------------------------------------------------------------------
\r
685 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
\r
687 // Give the value of cut set for the variable iVar and the pt bin iPtBin
\r
690 cout<<"Cuts not iniziaisez yet"<<endl;
\r
693 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
\r
695 //-------------------------------------------------------------------
\r
696 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
\r
698 // Get centrality percentile
\r
701 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
\r
702 if(mcArray) {fUseAOD049=kFALSE;}
\r
704 AliAODHeader *header=aodEvent->GetHeader();
\r
705 AliCentrality *centrality=header->GetCentralityP();
\r
706 Float_t cent=-999.;
\r
707 Bool_t isSelRun=kFALSE;
\r
708 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
\r
709 if(!centrality) return cent;
\r
711 if (estimator==kCentV0M){
\r
712 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
\r
714 Int_t quality = centrality->GetQuality();
\r
716 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
\r
718 Int_t runnum=aodEvent->GetRunNumber();
\r
719 for(Int_t ir=0;ir<5;ir++){
\r
720 if(runnum==selRun[ir]){
\r
725 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
\r
729 //temporary fix for AOD049 outliers
\r
730 if(fUseAOD049&¢>=0){
\r
732 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
\r
733 v0+=aodV0->GetMTotV0A();
\r
734 v0+=aodV0->GetMTotV0C();
\r
735 if(cent==0&&v0<19500)return -1;//filtering issue
\r
736 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
\r
737 Float_t val= 1.30552 + 0.147931 * v0;
\r
738 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};
\r
739 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
\r
743 if (estimator==kCentTRK) {
\r
744 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
\r
746 Int_t quality = centrality->GetQuality();
\r
748 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
\r
750 Int_t runnum=aodEvent->GetRunNumber();
\r
751 for(Int_t ir=0;ir<5;ir++){
\r
752 if(runnum==selRun[ir]){
\r
757 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
\r
762 if (estimator==kCentTKL){
\r
763 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
\r
765 Int_t quality = centrality->GetQuality();
\r
767 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
\r
769 Int_t runnum=aodEvent->GetRunNumber();
\r
770 for(Int_t ir=0;ir<5;ir++){
\r
771 if(runnum==selRun[ir]){
\r
776 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
\r
781 if (estimator==kCentCL1){
\r
782 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
\r
784 Int_t quality = centrality->GetQuality();
\r
786 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
\r
788 Int_t runnum=aodEvent->GetRunNumber();
\r
789 for(Int_t ir=0;ir<5;ir++){
\r
790 if(runnum==selRun[ir]){
\r
795 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
\r
800 AliWarning("Centrality estimator not valid");
\r
809 //-------------------------------------------------------------------
\r
810 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
\r
812 // Compare two cuts objects
\r
815 Bool_t areEqual=kTRUE;
\r
817 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
\r
819 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
\r
821 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
\r
823 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
\r
825 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
\r
827 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
\r
829 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
\r
831 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
\r
833 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
\r
835 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;}
\r
839 for(Int_t iv=0;iv<fnVars;iv++) {
\r
840 for(Int_t ib=0;ib<fnPtBins;ib++) {
\r
841 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
\r
842 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
\r
851 //---------------------------------------------------------------------------
\r
852 void AliRDHFCuts::MakeTable() const {
\r
854 // print cuts values in table format
\r
857 TString ptString = "pT range";
\r
858 if(fVarNames && fPtBinLimits && fCutsRD){
\r
859 TString firstLine(Form("* %-15s",ptString.Data()));
\r
860 for (Int_t ivar=0; ivar<fnVars; ivar++){
\r
861 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
\r
862 if (ivar == fnVars){
\r
866 Printf("%s",firstLine.Data());
\r
868 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
\r
870 if (ipt==fnPtBins-1){
\r
871 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
\r
874 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
\r
876 for (Int_t ivar=0; ivar<fnVars; ivar++){
\r
877 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
\r
879 Printf("%s",line.Data());
\r
887 //--------------------------------------------------------------------------
\r
888 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
\r
889 AliAODEvent *aod) const
\r
892 // Recalculate primary vertex without daughters
\r
896 AliError("Can not remove daughters from vertex without AOD event");
\r
900 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
\r
902 AliDebug(2,"Removal of daughter tracks failed");
\r
907 //set recalculed primary vertex
\r
908 d->SetOwnPrimaryVtx(recvtx);
\r
913 //--------------------------------------------------------------------------
\r
914 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
\r
917 // Recalculate primary vertex without daughters
\r
921 AliError("Can not get MC vertex without AOD event");
\r
926 AliAODMCHeader *mcHeader =
\r
927 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
\r
929 AliError("Can not get MC vertex without AODMCHeader event");
\r
933 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
\r
934 mcHeader->GetVertex(pos);
\r
935 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
\r
938 AliDebug(2,"Removal of daughter tracks failed");
\r
942 //set recalculed primary vertex
\r
943 d->SetOwnPrimaryVtx(recvtx);
\r
945 d->RecalculateImpPars(recvtx,aod);
\r
951 //--------------------------------------------------------------------------
\r
952 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
\r
954 AliAODVertex *origownvtx) const
\r
957 // Clean-up own primary vertex if needed
\r
960 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
\r
961 d->UnsetOwnPrimaryVtx();
\r
963 d->SetOwnPrimaryVtx(origownvtx);
\r
964 delete origownvtx; origownvtx=NULL;
\r
966 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
\r
969 delete origownvtx; origownvtx=NULL;
\r
974 //--------------------------------------------------------------------------
\r
975 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
\r
978 // Checks if this candidate is matched to MC signal
\r
981 if(!aod) return kFALSE;
\r
983 // get the MC array
\r
984 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
\r
986 if(!mcArray) return kFALSE;
\r
989 Int_t label = d->MatchToMC(pdg,mcArray);
\r
992 //printf("MATCH!\n");
\r