]> 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
Coverity
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCuts.cxx
CommitLineData
e0c4f114 1/**************************************************************************\r
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\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
15\r
16/* $Id$ */\r
17\r
18/////////////////////////////////////////////////////////////\r
19//\r
20// Base class for cuts on AOD reconstructed heavy-flavour decay\r
21//\r
22// Author: A.Dainese, andrea.dainese@pd.infn.it\r
23/////////////////////////////////////////////////////////////\r
24#include <Riostream.h>\r
25\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
31#include "AliLog.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
45\r
46ClassImp(AliRDHFCuts)\r
47\r
48//--------------------------------------------------------------------------\r
49AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) : \r
50AliAnalysisCuts(name,title),\r
51fMinVtxType(3),\r
52fMinVtxContr(1),\r
53fMaxVtxRedChi2(1e6),\r
54fMaxVtxZ(10.),\r
55fMinSPDMultiplicity(0),\r
56fTriggerMask(0),\r
57fTriggerClass("CINT1"),\r
58fTrackCuts(0),\r
59fnPtBins(1),\r
60fnPtBinLimits(1),\r
61fPtBinLimits(0),\r
62fnVars(1),\r
63fVarNames(0),\r
64fnVarsForOpt(0),\r
65fVarsForOpt(0),\r
66fGlobalIndex(1),\r
67fCutsRD(0),\r
68fIsUpperCut(0),\r
69fUsePID(kFALSE),\r
70fUseAOD049(kFALSE),\r
71fPidHF(0),\r
72fWhyRejection(0),\r
73fEvRejectionBits(0),\r
74fRemoveDaughtersFromPrimary(kFALSE),\r
75fUseMCVertex(kFALSE),\r
76fUsePhysicsSelection(kTRUE),\r
77fOptPileup(0),\r
78fMinContrPileup(3),\r
79fMinDzPileup(0.6),\r
80fUseCentrality(0),\r
81fMinCentrality(0.),\r
82fMaxCentrality(100.),\r
83fFixRefs(kFALSE),\r
84fIsSelectedCuts(0),\r
85fIsSelectedPID(0),\r
86fMinPtCand(-1.),\r
87fMaxPtCand(100000.),\r
88fKeepSignalMC(kFALSE),\r
89fIsCandTrackSPDFirst(kFALSE),\r
90fMaxPtCandTrackSPDFirst(0.)\r
91{\r
92 //\r
93 // Default Constructor\r
94 //\r
95}\r
96//--------------------------------------------------------------------------\r
97AliRDHFCuts::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
106 fTrackCuts(0),\r
107 fnPtBins(source.fnPtBins),\r
108 fnPtBinLimits(source.fnPtBinLimits),\r
109 fPtBinLimits(0),\r
110 fnVars(source.fnVars),\r
111 fVarNames(0),\r
112 fnVarsForOpt(source.fnVarsForOpt),\r
113 fVarsForOpt(0),\r
114 fGlobalIndex(source.fGlobalIndex),\r
115 fCutsRD(0),\r
116 fIsUpperCut(0),\r
117 fUsePID(source.fUsePID),\r
118 fUseAOD049(source.fUseAOD049),\r
119 fPidHF(0),\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
139{\r
140 //\r
141 // Copy constructor\r
142 //\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
150 PrintAll();\r
151\r
152}\r
153//--------------------------------------------------------------------------\r
154AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)\r
155{\r
156 //\r
157 // assignment operator\r
158 //\r
159 if(&source == this) return *this;\r
160\r
161 AliAnalysisCuts::operator=(source);\r
162\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
196\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
202 PrintAll();\r
203\r
204 return *this;\r
205}\r
206//--------------------------------------------------------------------------\r
207AliRDHFCuts::~AliRDHFCuts() {\r
208 // \r
209 // Default Destructor\r
210 //\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
214 if(fCutsRD) {\r
215 delete [] fCutsRD;\r
216 fCutsRD=0;\r
217 }\r
218 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}\r
219 if(fPidHF){ \r
220 delete fPidHF;\r
221 fPidHF=0;\r
222 }\r
223}\r
224//---------------------------------------------------------------------------\r
225Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {\r
226 //\r
227 // Centrality selection\r
228 //\r
229 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){ \r
230 AliWarning("Centrality estimator not valid"); \r
231 return 3; \r
232 }else{ \r
233 Float_t centvalue=GetCentrality((AliAODEvent*)event); \r
234 if (centvalue<-998.){//-999 if no centralityP\r
235 return 0; \r
236 }else{ \r
237 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){\r
238 return 2; \r
239 } \r
240 } \r
241 } \r
242 return 0;\r
243}\r
244//---------------------------------------------------------------------------\r
245Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {\r
246 //\r
247 // Event selection\r
248 // \r
249 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;\r
250\r
251 fWhyRejection=0;\r
252 fEvRejectionBits=0;\r
253 Bool_t accept=kTRUE;\r
254\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
259\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
269 }\r
270 // PbPb LHC10h\r
271 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);\r
272 // MC\r
273 if(isMC) fPidHF->SetMC(kTRUE);\r
274 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))\r
275 fPidHF->SetMClowenpp2011(kTRUE);\r
276 } \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
283 }\r
284 }\r
285\r
286\r
287 // trigger class\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
292 fWhyRejection=5;\r
293 fEvRejectionBits+=1<<kNotSelTrigger;\r
294 accept=kFALSE;\r
295 }\r
296 }\r
297\r
298 // TEMPORARY FIX FOR REFERENCES\r
299 // Fix references to daughter tracks\r
300 // if(fFixRefs) {\r
301 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();\r
302 // fixer->FixReferences((AliAODEvent*)event);\r
303 // delete fixer;\r
304 // }\r
305 //\r
306\r
307\r
308 // physics selection requirements\r
309 if(fUsePhysicsSelection){\r
310 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
311 if(!isSelected) {\r
312 if(accept) fWhyRejection=7;\r
313 fEvRejectionBits+=1<<kPhysicsSelection;\r
314 accept=kFALSE;\r
315 }\r
316 }\r
317\r
318 // vertex requirements\r
319 \r
320 const AliVVertex *vertex = event->GetPrimaryVertex();\r
321\r
322 if(!vertex){\r
323 accept=kFALSE;\r
324 fEvRejectionBits+=1<<kNoVertex;\r
325 }else{\r
326 TString title=vertex->GetTitle();\r
327 if(title.Contains("Z") && fMinVtxType>1){\r
328 accept=kFALSE;\r
329 fEvRejectionBits+=1<<kNoVertex;\r
330 }\r
331 else if(title.Contains("3D") && fMinVtxType>2){\r
332 accept=kFALSE;\r
333 fEvRejectionBits+=1<<kNoVertex;\r
334 }\r
335 if(vertex->GetNContributors()<fMinVtxContr){\r
336 accept=kFALSE;\r
337 fEvRejectionBits+=1<<kTooFewVtxContrib;\r
338 }\r
339 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {\r
340 fEvRejectionBits+=1<<kZVtxOutFid;\r
341 if(accept) fWhyRejection=6;\r
342 accept=kFALSE;\r
343 } \r
344 }\r
345\r
346\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
354 accept=kFALSE;\r
355 }\r
356 }\r
357\r
358 // centrality selection\r
359 if (fUseCentrality!=kCentOff) { \r
360 Int_t rejection=IsEventSelectedInCentrality(event); \r
361 if(rejection>1){ \r
362 if(accept) fWhyRejection=rejection; \r
363 fEvRejectionBits+=1<<kOutsideCentrality;\r
364 accept=kFALSE;\r
365 }\r
366 }\r
367\r
368 return accept;\r
369}\r
370//---------------------------------------------------------------------------\r
371Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {\r
372 //\r
373 // Daughter tracks selection\r
374 // \r
375 if(!fTrackCuts) return kTRUE;\r
376 \r
377 Int_t ndaughters = d->GetNDaughters();\r
378 AliAODVertex *vAOD = d->GetPrimaryVtx();\r
379 Double_t pos[3],cov[6];\r
380 vAOD->GetXYZ(pos);\r
381 vAOD->GetCovarianceMatrix(cov);\r
382 const AliESDVertex vESD(pos,cov,100.,100);\r
383 \r
384 Bool_t retval=kTRUE;\r
385 \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
391\r
392 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)\r
393 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }\r
394\r
395 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;\r
396 }\r
397 \r
398 return retval;\r
399}\r
400//---------------------------------------------------------------------------\r
401Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {\r
402 //\r
403 // Convert to ESDtrack, relate to vertex and check cuts\r
404 //\r
405 if(!cuts) return kTRUE;\r
406\r
407 Bool_t retval=kTRUE;\r
408\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
417\r
418 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;\r
419 \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
424 }\r
425 return retval; \r
426}\r
427//---------------------------------------------------------------------------\r
428void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {\r
429 // Set the pt bins\r
430\r
431 if(fPtBinLimits) {\r
432 delete [] fPtBinLimits;\r
433 fPtBinLimits = NULL;\r
434 printf("Changing the pt bins\n");\r
435 }\r
436\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
440 }\r
441\r
442 fnPtBinLimits = nPtBinLimits;\r
443 SetGlobalIndex();\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
447\r
448 return;\r
449}\r
450//---------------------------------------------------------------------------\r
451void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){\r
452 // Set the variable names\r
453\r
454 if(fVarNames) {\r
455 delete [] fVarNames;\r
456 fVarNames = NULL;\r
457 //printf("Changing the variable names\n");\r
458 }\r
459 if(nVars!=fnVars){\r
460 printf("Wrong number of variables: it has to be %d\n",fnVars);\r
461 return;\r
462 }\r
463 //fnVars=nVars;\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
469 }\r
470\r
471 return;\r
472}\r
473//---------------------------------------------------------------------------\r
474void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {\r
475 // Set the variables to be used for cuts optimization\r
476\r
477 if(fVarsForOpt) {\r
478 delete [] fVarsForOpt;\r
479 fVarsForOpt = NULL;\r
480 //printf("Changing the variables for cut optimization\n");\r
481 }\r
482 \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
485 return;\r
486 } \r
487 \r
488 fnVarsForOpt = 0;\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
493 }\r
494\r
495 return;\r
496}\r
497\r
498//---------------------------------------------------------------------------\r
499void AliRDHFCuts::SetUseCentrality(Int_t flag) {\r
500 //\r
501 // set centrality estimator \r
502 //\r
503 fUseCentrality=flag;\r
504 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");\r
505 \r
506 return;\r
507}\r
508\r
509\r
510//---------------------------------------------------------------------------\r
511void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {\r
512 //\r
513 // store the cuts\r
514 //\r
515 if(nVars!=fnVars) {\r
516 printf("Wrong number of variables: it has to be %d\n",fnVars);\r
517 AliFatal("exiting");\r
518 } \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
522 } \r
523\r
524 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];\r
525 \r
526\r
527 for(Int_t iv=0; iv<fnVars; iv++) {\r
528\r
529 for(Int_t ib=0; ib<fnPtBins; ib++) {\r
530\r
531 //check\r
532 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {\r
533 cout<<"Overflow, exit..."<<endl;\r
534 return;\r
535 }\r
536\r
537 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];\r
538\r
539 }\r
540 }\r
541 return;\r
542}\r
543//---------------------------------------------------------------------------\r
544void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){\r
545 //\r
546 // store the cuts\r
547 //\r
548 if(glIndex != fGlobalIndex){\r
549 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;\r
550 AliFatal("exiting");\r
551 }\r
552 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];\r
553\r
554 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){\r
555 fCutsRD[iGl] = cutsRDGlob[iGl];\r
556 }\r
557 return;\r
558}\r
559//---------------------------------------------------------------------------\r
560void AliRDHFCuts::PrintAll() const {\r
561 //\r
562 // print all cuts values\r
563 // \r
564\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
581 }\r
582 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);\r
583\r
584 if(fVarNames){\r
585 cout<<"Array of variables"<<endl;\r
586 for(Int_t iv=0;iv<fnVars;iv++){\r
587 cout<<fVarNames[iv]<<"\t";\r
588 }\r
589 cout<<endl;\r
590 }\r
591 if(fVarsForOpt){\r
592 cout<<"Array of optimization"<<endl;\r
593 for(Int_t iv=0;iv<fnVars;iv++){\r
594 cout<<fVarsForOpt[iv]<<"\t";\r
595 }\r
596 cout<<endl;\r
597 }\r
598 if(fIsUpperCut){\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
602 }\r
603 cout<<endl;\r
604 }\r
605 if(fPtBinLimits){\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
609 }\r
610 cout<<endl;\r
611 }\r
612 if(fCutsRD){\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
617 } \r
618 cout<<endl;\r
619 }\r
620 cout<<endl;\r
621 }\r
622 return;\r
623}\r
624//---------------------------------------------------------------------------\r
625void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{\r
626 //\r
627 // get the cuts\r
628 //\r
629\r
630 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;\r
631\r
632\r
633 Int_t iv,ib;\r
634 if(!cutsRD) {\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
639 }\r
640 }\r
641 \r
642 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {\r
643 GetVarPtIndex(iGlobal,iv,ib);\r
644 cutsRD[iv][ib] = fCutsRD[iGlobal];\r
645 }\r
646\r
647 return;\r
648}\r
649\r
650//---------------------------------------------------------------------------\r
651Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{\r
652 //\r
653 // give the global index from variable and pt bin\r
654 //\r
655 return iPtBin*fnVars+iVar;\r
656}\r
657\r
658//---------------------------------------------------------------------------\r
659void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {\r
660 //\r
661 //give the index of the variable and of the pt bin from the global index\r
662 //\r
663 iPtBin=(Int_t)iGlob/fnVars;\r
664 iVar=iGlob%fnVars;\r
665\r
666 return;\r
667}\r
668\r
669//---------------------------------------------------------------------------\r
670Int_t AliRDHFCuts::PtBin(Double_t pt) const {\r
671 //\r
672 //give the pt bin where the pt lies.\r
673 //\r
674 Int_t ptbin=-1;\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
678 ptbin=i;\r
679 break;\r
680 }\r
681 }\r
682 return ptbin;\r
683}\r
684//-------------------------------------------------------------------\r
685Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {\r
686 // \r
687 // Give the value of cut set for the variable iVar and the pt bin iPtBin\r
688 //\r
689 if(!fCutsRD){\r
690 cout<<"Cuts not iniziaisez yet"<<endl;\r
691 return 0;\r
692 }\r
693 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];\r
694}\r
695//-------------------------------------------------------------------\r
696Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {\r
697 //\r
698 // Get centrality percentile\r
699 //\r
700\r
701 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
702 if(mcArray) {fUseAOD049=kFALSE;}\r
703\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
710 else{\r
711 if (estimator==kCentV0M){\r
712 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));\r
713 if(cent<0){\r
714 Int_t quality = centrality->GetQuality();\r
715 if(quality<=1){\r
716 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
717 }else{\r
718 Int_t runnum=aodEvent->GetRunNumber();\r
719 for(Int_t ir=0;ir<5;ir++){\r
720 if(runnum==selRun[ir]){\r
721 isSelRun=kTRUE;\r
722 break;\r
723 }\r
724 }\r
725 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");\r
726 }\r
727 }\r
728\r
729 //temporary fix for AOD049 outliers\r
730 if(fUseAOD049&&cent>=0){\r
731 Float_t v0=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
740 }\r
741 }\r
742 else {\r
743 if (estimator==kCentTRK) {\r
744 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));\r
745 if(cent<0){\r
746 Int_t quality = centrality->GetQuality();\r
747 if(quality<=1){\r
748 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
749 }else{\r
750 Int_t runnum=aodEvent->GetRunNumber();\r
751 for(Int_t ir=0;ir<5;ir++){\r
752 if(runnum==selRun[ir]){\r
753 isSelRun=kTRUE;\r
754 break;\r
755 }\r
756 }\r
757 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");\r
758 }\r
759 }\r
760 }\r
761 else{\r
762 if (estimator==kCentTKL){\r
763 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));\r
764 if(cent<0){\r
765 Int_t quality = centrality->GetQuality();\r
766 if(quality<=1){\r
767 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
768 }else{\r
769 Int_t runnum=aodEvent->GetRunNumber();\r
770 for(Int_t ir=0;ir<5;ir++){\r
771 if(runnum==selRun[ir]){\r
772 isSelRun=kTRUE;\r
773 break;\r
774 }\r
775 }\r
776 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");\r
777 } \r
778 }\r
779 }\r
780 else{\r
781 if (estimator==kCentCL1){\r
782 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));\r
783 if(cent<0){\r
784 Int_t quality = centrality->GetQuality();\r
785 if(quality<=1){\r
786 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
787 }else{\r
788 Int_t runnum=aodEvent->GetRunNumber();\r
789 for(Int_t ir=0;ir<5;ir++){\r
790 if(runnum==selRun[ir]){\r
791 isSelRun=kTRUE;\r
792 break;\r
793 }\r
794 }\r
795 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");\r
796 }\r
797 }\r
798 }\r
799 else {\r
800 AliWarning("Centrality estimator not valid");\r
801 \r
802 }\r
803 }\r
804 }\r
805 } \r
806 }\r
807 return cent;\r
808}\r
809//-------------------------------------------------------------------\r
810Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {\r
811 //\r
812 // Compare two cuts objects\r
813 //\r
814\r
815 Bool_t areEqual=kTRUE;\r
816\r
817 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}\r
818\r
819 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}\r
820\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
822\r
823 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}\r
824\r
825 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}\r
826\r
827 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}\r
828 if(fTrackCuts){\r
829 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}\r
830\r
831 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}\r
832\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
834\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
836 }\r
837\r
838 if(fCutsRD) {\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
843 areEqual=kFALSE;\r
844 }\r
845 }\r
846 }\r
847 }\r
848\r
849 return areEqual;\r
850}\r
851//---------------------------------------------------------------------------\r
852void AliRDHFCuts::MakeTable() const {\r
853 //\r
854 // print cuts values in table format\r
855 // \r
856\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
863 firstLine+="*\n";\r
864 }\r
865 }\r
866 Printf("%s",firstLine.Data());\r
867 \r
868 for (Int_t ipt=0; ipt<fnPtBins; ipt++){\r
869 TString line;\r
870 if (ipt==fnPtBins-1){\r
871 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);\r
872 }\r
873 else{\r
874 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);\r
875 }\r
876 for (Int_t ivar=0; ivar<fnVars; ivar++){\r
877 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);\r
878 }\r
879 Printf("%s",line.Data());\r
880 }\r
881\r
882 }\r
883\r
884\r
885 return;\r
886}\r
887//--------------------------------------------------------------------------\r
888Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
889 AliAODEvent *aod) const\r
890{\r
891 //\r
892 // Recalculate primary vertex without daughters\r
893 //\r
894\r
895 if(!aod) {\r
896 AliError("Can not remove daughters from vertex without AOD event");\r
897 return 0;\r
898 } \r
899\r
900 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);\r
901 if(!recvtx){\r
902 AliDebug(2,"Removal of daughter tracks failed");\r
903 return kFALSE;\r
904 }\r
905\r
906\r
907 //set recalculed primary vertex\r
908 d->SetOwnPrimaryVtx(recvtx);\r
909 delete recvtx;\r
910\r
911 return kTRUE;\r
912}\r
913//--------------------------------------------------------------------------\r
914Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const\r
915{\r
916 //\r
917 // Recalculate primary vertex without daughters\r
918 //\r
919\r
920 if(!aod) {\r
921 AliError("Can not get MC vertex without AOD event");\r
922 return kFALSE;\r
923 } \r
924\r
925 // load MC header\r
926 AliAODMCHeader *mcHeader = \r
927 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
928 if(!mcHeader) {\r
929 AliError("Can not get MC vertex without AODMCHeader event");\r
930 return kFALSE;\r
931 }\r
932 Double_t pos[3];\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
936\r
937 if(!recvtx){\r
938 AliDebug(2,"Removal of daughter tracks failed");\r
939 return kFALSE;\r
940 }\r
941\r
942 //set recalculed primary vertex\r
943 d->SetOwnPrimaryVtx(recvtx);\r
944\r
945 d->RecalculateImpPars(recvtx,aod);\r
946\r
947 delete recvtx;\r
948\r
949 return kTRUE;\r
950}\r
951//--------------------------------------------------------------------------\r
952void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,\r
953 AliAODEvent *aod,\r
954 AliAODVertex *origownvtx) const\r
955{\r
956 //\r
957 // Clean-up own primary vertex if needed\r
958 //\r
959\r
960 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {\r
961 d->UnsetOwnPrimaryVtx();\r
962 if(origownvtx) {\r
963 d->SetOwnPrimaryVtx(origownvtx);\r
964 delete origownvtx; origownvtx=NULL;\r
965 }\r
966 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);\r
967 } else {\r
968 if(origownvtx) {\r
969 delete origownvtx; origownvtx=NULL;\r
970 }\r
971 }\r
972 return;\r
973}\r
974//--------------------------------------------------------------------------\r
975Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const \r
976{\r
977 //\r
978 // Checks if this candidate is matched to MC signal\r
979 // \r
980\r
981 if(!aod) return kFALSE;\r
982\r
983 // get the MC array\r
984 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
985\r
986 if(!mcArray) return kFALSE;\r
987\r
988 // try to match \r
989 Int_t label = d->MatchToMC(pdg,mcArray);\r
990 \r
991 if(label>=0) {\r
992 //printf("MATCH!\n");\r
993 return kTRUE;\r
994 }\r
995\r
996 return kFALSE;\r
997}\r
998\r
999\r