]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliRDHFCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o updates for AOD MC (Julian)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCuts.cxx
CommitLineData
0c22e2ac 1/**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18/////////////////////////////////////////////////////////////
19//
20// Base class for cuts on AOD reconstructed heavy-flavour decay
21//
22// Author: A.Dainese, andrea.dainese@pd.infn.it
23/////////////////////////////////////////////////////////////
24#include <Riostream.h>
25
26#include "AliVEvent.h"
27#include "AliESDEvent.h"
28#include "AliAODEvent.h"
29#include "AliVVertex.h"
30#include "AliESDVertex.h"
31#include "AliLog.h"
32#include "AliAODVertex.h"
33#include "AliESDtrack.h"
34#include "AliAODTrack.h"
35#include "AliESDtrackCuts.h"
36#include "AliCentrality.h"
37#include "AliAODRecoDecayHF.h"
38#include "AliAnalysisVertexingHF.h"
39#include "AliAODMCHeader.h"
40#include "AliAODMCParticle.h"
4bdc0315 41#include "AliVertexerTracks.h"
0c22e2ac 42#include "AliRDHFCuts.h"
43#include "AliAnalysisManager.h"
44#include "AliInputEventHandler.h"
45#include "AliPIDResponse.h"
46
47ClassImp(AliRDHFCuts)
48
49//--------------------------------------------------------------------------
50AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
51AliAnalysisCuts(name,title),
52fMinVtxType(3),
53fMinVtxContr(1),
54fMaxVtxRedChi2(1e6),
55fMaxVtxZ(10.),
56fMinSPDMultiplicity(0),
57fTriggerMask(AliVEvent::kAnyINT),
58fUseOnlyOneTrigger(kFALSE),
0c22e2ac 59fTrackCuts(0),
60fnPtBins(1),
61fnPtBinLimits(1),
62fPtBinLimits(0),
63fnVars(1),
64fVarNames(0),
65fnVarsForOpt(0),
66fVarsForOpt(0),
67fGlobalIndex(1),
68fCutsRD(0),
69fIsUpperCut(0),
70fUsePID(kFALSE),
71fUseAOD049(kFALSE),
72fPidHF(0),
73fWhyRejection(0),
74fEvRejectionBits(0),
75fRemoveDaughtersFromPrimary(kFALSE),
4bdc0315 76fRecomputePrimVertex(kFALSE),
0c22e2ac 77fUseMCVertex(kFALSE),
78fUsePhysicsSelection(kTRUE),
79fOptPileup(0),
80fMinContrPileup(3),
81fMinDzPileup(0.6),
82fUseCentrality(0),
83fMinCentrality(0.),
84fMaxCentrality(100.),
85fFixRefs(kFALSE),
86fIsSelectedCuts(0),
87fIsSelectedPID(0),
88fMinPtCand(-1.),
89fMaxPtCand(100000.),
90fKeepSignalMC(kFALSE),
91fIsCandTrackSPDFirst(kFALSE),
851d1cae 92fMaxPtCandTrackSPDFirst(0.),
93fApplySPDDeadPbPb2011(kFALSE),
fdd5a897 94fRemoveTrackletOutliers(kFALSE),
95fKinkReject(kFALSE)
0c22e2ac 96{
97 //
98 // Default Constructor
99 //
68c15f66 100 fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
0c22e2ac 101}
102//--------------------------------------------------------------------------
103AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
104 AliAnalysisCuts(source),
105 fMinVtxType(source.fMinVtxType),
106 fMinVtxContr(source.fMinVtxContr),
107 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
108 fMaxVtxZ(source.fMaxVtxZ),
109 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
110 fTriggerMask(source.fTriggerMask),
111 fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
68c15f66 112 fTriggerClass(),
0c22e2ac 113 fTrackCuts(0),
114 fnPtBins(source.fnPtBins),
115 fnPtBinLimits(source.fnPtBinLimits),
116 fPtBinLimits(0),
117 fnVars(source.fnVars),
118 fVarNames(0),
119 fnVarsForOpt(source.fnVarsForOpt),
120 fVarsForOpt(0),
121 fGlobalIndex(source.fGlobalIndex),
122 fCutsRD(0),
123 fIsUpperCut(0),
124 fUsePID(source.fUsePID),
125 fUseAOD049(source.fUseAOD049),
126 fPidHF(0),
127 fWhyRejection(source.fWhyRejection),
128 fEvRejectionBits(source.fEvRejectionBits),
129 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
4bdc0315 130 fRecomputePrimVertex(source.fRecomputePrimVertex),
0c22e2ac 131 fUseMCVertex(source.fUseMCVertex),
132 fUsePhysicsSelection(source.fUsePhysicsSelection),
133 fOptPileup(source.fOptPileup),
134 fMinContrPileup(source.fMinContrPileup),
135 fMinDzPileup(source.fMinDzPileup),
136 fUseCentrality(source.fUseCentrality),
137 fMinCentrality(source.fMinCentrality),
138 fMaxCentrality(source.fMaxCentrality),
139 fFixRefs(source.fFixRefs),
140 fIsSelectedCuts(source.fIsSelectedCuts),
141 fIsSelectedPID(source.fIsSelectedPID),
142 fMinPtCand(source.fMinPtCand),
143 fMaxPtCand(source.fMaxPtCand),
144 fKeepSignalMC(source.fKeepSignalMC),
145 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
851d1cae 146 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
147 fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
fdd5a897 148 fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
149 fKinkReject(source.fKinkReject)
0c22e2ac 150{
151 //
152 // Copy constructor
153 //
154 cout<<"Copy constructor"<<endl;
68c15f66 155 fTriggerClass[0] = source.fTriggerClass[0];
156 fTriggerClass[1] = source.fTriggerClass[1];
0c22e2ac 157 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
158 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
159 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
160 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
161 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
162 if(source.fPidHF) SetPidHF(source.fPidHF);
163 PrintAll();
164
165}
166//--------------------------------------------------------------------------
167AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
168{
169 //
170 // assignment operator
171 //
172 if(&source == this) return *this;
173
174 AliAnalysisCuts::operator=(source);
175
176 fMinVtxType=source.fMinVtxType;
177 fMinVtxContr=source.fMinVtxContr;
178 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
179 fMaxVtxZ=source.fMaxVtxZ;
180 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
181 fTriggerMask=source.fTriggerMask;
182 fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
68c15f66 183 fTriggerClass[0]=source.fTriggerClass[0];
184 fTriggerClass[1]=source.fTriggerClass[1];
0c22e2ac 185 fnPtBins=source.fnPtBins;
186 fnPtBinLimits=source.fnPtBinLimits;
187 fnVars=source.fnVars;
188 fGlobalIndex=source.fGlobalIndex;
189 fnVarsForOpt=source.fnVarsForOpt;
190 fUsePID=source.fUsePID;
191 fUseAOD049=source.fUseAOD049;
192 if(fPidHF) delete fPidHF;
193 fPidHF=new AliAODPidHF(*(source.GetPidHF()));
194 fWhyRejection=source.fWhyRejection;
195 fEvRejectionBits=source.fEvRejectionBits;
196 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
4bdc0315 197 fRecomputePrimVertex=source.fRecomputePrimVertex;
0c22e2ac 198 fUseMCVertex=source.fUseMCVertex;
199 fUsePhysicsSelection=source.fUsePhysicsSelection;
200 fOptPileup=source.fOptPileup;
201 fMinContrPileup=source.fMinContrPileup;
202 fMinDzPileup=source.fMinDzPileup;
203 fUseCentrality=source.fUseCentrality;
204 fMinCentrality=source.fMinCentrality;
205 fMaxCentrality=source.fMaxCentrality;
206 fFixRefs=source.fFixRefs;
207 fIsSelectedCuts=source.fIsSelectedCuts;
208 fIsSelectedPID=source.fIsSelectedPID;
209 fMinPtCand=source.fMinPtCand;
210 fMaxPtCand=source.fMaxPtCand;
211 fKeepSignalMC=source.fKeepSignalMC;
212 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
213 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
851d1cae 214 fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
215 fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
fdd5a897 216 fKinkReject=source.fKinkReject;
0c22e2ac 217
218 if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
219 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
220 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
221 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
222 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
223 PrintAll();
224
225 return *this;
226}
227//--------------------------------------------------------------------------
228AliRDHFCuts::~AliRDHFCuts() {
229 //
230 // Default Destructor
231 //
232 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
233 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
234 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
235 if(fCutsRD) {
236 delete [] fCutsRD;
237 fCutsRD=0;
238 }
239 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
240 if(fPidHF){
241 delete fPidHF;
242 fPidHF=0;
243 }
244}
245//---------------------------------------------------------------------------
246Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
247 //
248 // Centrality selection
249 //
250 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
251 AliWarning("Centrality estimator not valid");
252 return 3;
253 }else{
254 Float_t centvalue=GetCentrality((AliAODEvent*)event);
255 if (centvalue<-998.){//-999 if no centralityP
256 return 0;
257 }else{
258 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
259 return 2;
260 }
261 }
262 }
263 return 0;
264}
265//---------------------------------------------------------------------------
266Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
267 //
268 // Event selection
269 //
270 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
271
272 fWhyRejection=0;
273 fEvRejectionBits=0;
274 Bool_t accept=kTRUE;
275
4bdc0315 276 if(fRecomputePrimVertex){
277 Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
278 if(!vertOK){
279 fWhyRejection=6;
280 return kFALSE;
281 }
282 }
283
0c22e2ac 284 // check if it's MC
285 Bool_t isMC=kFALSE;
286 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
287 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
288
289 // settings for the TPC dE/dx BB parameterization
347d8758 290 if(fPidHF){
291 if(fPidHF->GetPidResponse()==0x0){
292 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
293 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
294 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
295 fPidHF->SetPidResponse(pidResp);
296 }
297 if(fPidHF->GetOldPid()) {
298 // pp, from LHC10d onwards
299 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
300 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
301 // pp, 2011 low energy run
302 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
303 fPidHF->SetppLowEn2011(kTRUE);
304 fPidHF->SetOnePad(kFALSE);
305 }
306 // PbPb LHC10h
307 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
308 // MC
309 if(isMC) fPidHF->SetMC(kTRUE);
310 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
311 fPidHF->SetMClowenpp2011(kTRUE);
312 fPidHF->SetBetheBloch();
313 }else{
314 // check that AliPIDResponse object was properly set in case of using OADB
315 if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
0c22e2ac 316 }
0c22e2ac 317 }
318
319
320 // trigger class
321 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
322 // don't do for MC and for PbPb 2010 data
323 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
68c15f66 324 if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
325 (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
0c22e2ac 326 fWhyRejection=5;
327 fEvRejectionBits+=1<<kNotSelTrigger;
328 accept=kFALSE;
329 }
330 }
331
332 // TEMPORARY FIX FOR REFERENCES
333 // Fix references to daughter tracks
334 // if(fFixRefs) {
335 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
336 // fixer->FixReferences((AliAODEvent*)event);
337 // delete fixer;
338 // }
339 //
340
341
342 // physics selection requirements
343 if(fUsePhysicsSelection){
344 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
345 if(!isSelected) {
346 if(accept) fWhyRejection=7;
347 fEvRejectionBits+=1<<kPhysicsSelection;
348 accept=kFALSE;
349 }else{
350 if(fUseOnlyOneTrigger){
351 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
352 if(accept) fWhyRejection=7;
353 fEvRejectionBits+=1<<kPhysicsSelection;
354 accept=kFALSE;
355 }
356 }
357 }
358 }
359
360 // vertex requirements
361
362 const AliVVertex *vertex = event->GetPrimaryVertex();
363
364 if(!vertex){
365 accept=kFALSE;
366 fEvRejectionBits+=1<<kNoVertex;
367 }else{
368 TString title=vertex->GetTitle();
369 if(title.Contains("Z") && fMinVtxType>1){
370 accept=kFALSE;
371 fEvRejectionBits+=1<<kNoVertex;
372 }
373 else if(title.Contains("3D") && fMinVtxType>2){
374 accept=kFALSE;
375 fEvRejectionBits+=1<<kNoVertex;
376 }
377 if(vertex->GetNContributors()<fMinVtxContr){
378 accept=kFALSE;
379 fEvRejectionBits+=1<<kTooFewVtxContrib;
380 }
381 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
382 fEvRejectionBits+=1<<kZVtxOutFid;
383 if(accept) fWhyRejection=6;
384 accept=kFALSE;
385 }
386 }
387
388
389 // pile-up rejection
390 if(fOptPileup==kRejectPileupEvent){
391 Int_t cutc=(Int_t)fMinContrPileup;
392 Double_t cutz=(Double_t)fMinDzPileup;
393 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
394 if(accept) fWhyRejection=1;
395 fEvRejectionBits+=1<<kPileupSPD;
396 accept=kFALSE;
397 }
398 }
399
400 // centrality selection
401 if (fUseCentrality!=kCentOff) {
402 Int_t rejection=IsEventSelectedInCentrality(event);
403 if(rejection>1){
404 if(accept) fWhyRejection=rejection;
405 fEvRejectionBits+=1<<kOutsideCentrality;
406 accept=kFALSE;
407 }
408 }
409
851d1cae 410 // PbPb2011 outliers in tracklets vs. VZERO
411 if(!isMC && event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
412 if(fRemoveTrackletOutliers){
413 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
414 Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
415 Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
416 if(v0cent<cutval){
417 if(accept) fWhyRejection=2;
418 fEvRejectionBits+=1<<kOutsideCentrality;
419 accept=kFALSE;
420 }
421 }
422 }
423
0c22e2ac 424 return accept;
425}
426//---------------------------------------------------------------------------
427Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
428 //
429 // Daughter tracks selection
430 //
431 if(!fTrackCuts) return kTRUE;
432
433 Int_t ndaughters = d->GetNDaughters();
434 AliAODVertex *vAOD = d->GetPrimaryVtx();
435 Double_t pos[3],cov[6];
436 vAOD->GetXYZ(pos);
437 vAOD->GetCovarianceMatrix(cov);
438 const AliESDVertex vESD(pos,cov,100.,100);
439
440 Bool_t retval=kTRUE;
441
442 for(Int_t idg=0; idg<ndaughters; idg++) {
443 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
444 if(!dgTrack) {retval = kFALSE; continue;}
445 //printf("charge %d\n",dgTrack->Charge());
446 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
447
448 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
449 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
450
451 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
452 }
453
454 return retval;
455}
456//---------------------------------------------------------------------------
457Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
458 //
459 // Convert to ESDtrack, relate to vertex and check cuts
460 //
461 if(!cuts) return kTRUE;
462
463 if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
464
465 Bool_t retval=kTRUE;
466
467 // convert to ESD track here
468 AliESDtrack esdTrack(track);
469 // set the TPC cluster info
470 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
471 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
472 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
473 // needed to calculate the impact parameters
474 esdTrack.RelateToVertex(primary,0.,3.);
475
476 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
fdd5a897 477
478
479 if(fKinkReject){
480 AliAODVertex *maybeKink=track->GetProdVertex();
481 if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
482 }
0c22e2ac 483
484 if(fOptPileup==kRejectTracksFromPileupVertex){
485 // to be implemented
486 // we need either to have here the AOD Event,
487 // or to have the pileup vertex object
488 }
851d1cae 489
490 if(fApplySPDDeadPbPb2011){
491 Bool_t deadSPDLay1PbPb2011[20][4]={
492 {kTRUE,kTRUE,kTRUE,kTRUE},
493 {kTRUE,kTRUE,kTRUE,kTRUE},
494 {kTRUE,kTRUE,kTRUE,kTRUE},
495 {kTRUE,kTRUE,kTRUE,kTRUE},
496 {kTRUE,kTRUE,kTRUE,kTRUE},
497 {kFALSE,kFALSE,kTRUE,kTRUE},
498 {kTRUE,kTRUE,kFALSE,kFALSE},
499 {kTRUE,kTRUE,kTRUE,kTRUE},
500 {kFALSE,kFALSE,kTRUE,kTRUE},
501 {kTRUE,kTRUE,kTRUE,kTRUE},
502 {kTRUE,kTRUE,kFALSE,kFALSE},
503 {kTRUE,kTRUE,kTRUE,kTRUE},
504 {kFALSE,kFALSE,kFALSE,kFALSE},
505 {kFALSE,kFALSE,kTRUE,kTRUE},
506 {kFALSE,kFALSE,kFALSE,kFALSE},
507 {kFALSE,kFALSE,kFALSE,kFALSE},
508 {kTRUE,kTRUE,kTRUE,kTRUE},
509 {kTRUE,kTRUE,kFALSE,kFALSE},
510 {kFALSE,kFALSE,kFALSE,kFALSE},
511 {kFALSE,kFALSE,kFALSE,kFALSE}
512 };
513 Bool_t deadSPDLay2PbPb2011[40][4]={
514 {kTRUE,kTRUE,kTRUE,kTRUE},
515 {kTRUE,kTRUE,kTRUE,kTRUE},
516 {kTRUE,kTRUE,kTRUE,kTRUE},
517 {kTRUE,kTRUE,kTRUE,kTRUE},
518 {kTRUE,kTRUE,kTRUE,kTRUE},
519 {kTRUE,kTRUE,kTRUE,kTRUE},
520 {kTRUE,kTRUE,kTRUE,kTRUE},
521 {kTRUE,kTRUE,kTRUE,kTRUE},
522 {kTRUE,kTRUE,kTRUE,kTRUE},
523 {kTRUE,kTRUE,kTRUE,kTRUE},
524 {kTRUE,kTRUE,kTRUE,kTRUE},
525 {kTRUE,kTRUE,kTRUE,kTRUE},
526 {kFALSE,kFALSE,kFALSE,kFALSE},
527 {kFALSE,kFALSE,kTRUE,kTRUE},
528 {kTRUE,kTRUE,kTRUE,kTRUE},
529 {kTRUE,kTRUE,kTRUE,kTRUE},
530 {kTRUE,kTRUE,kFALSE,kFALSE},
531 {kTRUE,kTRUE,kTRUE,kTRUE},
532 {kTRUE,kTRUE,kTRUE,kTRUE},
533 {kTRUE,kTRUE,kTRUE,kTRUE},
534 {kFALSE,kFALSE,kFALSE,kFALSE},
535 {kFALSE,kFALSE,kFALSE,kFALSE},
536 {kTRUE,kTRUE,kTRUE,kTRUE},
537 {kTRUE,kTRUE,kTRUE,kTRUE},
538 {kFALSE,kFALSE,kFALSE,kFALSE},
539 {kFALSE,kFALSE,kFALSE,kFALSE},
540 {kTRUE,kTRUE,kTRUE,kTRUE},
541 {kTRUE,kTRUE,kTRUE,kTRUE},
542 {kFALSE,kFALSE,kFALSE,kFALSE},
543 {kFALSE,kFALSE,kFALSE,kFALSE},
544 {kFALSE,kFALSE,kFALSE,kFALSE},
545 {kFALSE,kFALSE,kFALSE,kFALSE},
546 {kTRUE,kTRUE,kTRUE,kTRUE},
547 {kTRUE,kTRUE,kTRUE,kTRUE},
548 {kTRUE,kTRUE,kFALSE,kFALSE},
549 {kTRUE,kTRUE,kTRUE,kTRUE},
550 {kFALSE,kFALSE,kFALSE,kFALSE},
551 {kFALSE,kFALSE,kFALSE,kFALSE},
552 {kFALSE,kFALSE,kFALSE,kFALSE},
553 {kFALSE,kFALSE,kFALSE,kFALSE}
554 };
555 Double_t xyz1[3],xyz2[3];
556 esdTrack.GetXYZAt(3.9,0.,xyz1);
557 esdTrack.GetXYZAt(7.6,0.,xyz2);
558 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
559 if(phi1<0) phi1+=2*TMath::Pi();
a0dabf3d 560 Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
851d1cae 561 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
562 if(phi2<0) phi2+=2*TMath::Pi();
a0dabf3d 563 Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
564 Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
565 Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
851d1cae 566 Bool_t lay1ok=kFALSE;
567 if(mod1>=0 && mod1<4 && lad1<20){
568 lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
569 }
570 Bool_t lay2ok=kFALSE;
571 if(mod2>=0 && mod2<4 && lad2<40){
572 lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
573 }
574 if(!lay1ok && !lay2ok) retval=kFALSE;
575 }
576
0c22e2ac 577 return retval;
578}
579//---------------------------------------------------------------------------
580void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
581 // Set the pt bins
582
583 if(fPtBinLimits) {
584 delete [] fPtBinLimits;
585 fPtBinLimits = NULL;
586 printf("Changing the pt bins\n");
587 }
588
589 if(nPtBinLimits != fnPtBins+1){
590 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
591 SetNPtBins(nPtBinLimits-1);
592 }
593
594 fnPtBinLimits = nPtBinLimits;
595 SetGlobalIndex();
596 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
597 fPtBinLimits = new Float_t[fnPtBinLimits];
598 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
599
600 return;
601}
602//---------------------------------------------------------------------------
603void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
604 // Set the variable names
605
606 if(fVarNames) {
607 delete [] fVarNames;
608 fVarNames = NULL;
609 //printf("Changing the variable names\n");
610 }
611 if(nVars!=fnVars){
612 printf("Wrong number of variables: it has to be %d\n",fnVars);
613 return;
614 }
615 //fnVars=nVars;
616 fVarNames = new TString[nVars];
617 fIsUpperCut = new Bool_t[nVars];
618 for(Int_t iv=0; iv<nVars; iv++) {
619 fVarNames[iv] = varNames[iv];
620 fIsUpperCut[iv] = isUpperCut[iv];
621 }
622
623 return;
624}
625//---------------------------------------------------------------------------
626void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
627 // Set the variables to be used for cuts optimization
628
629 if(fVarsForOpt) {
630 delete [] fVarsForOpt;
631 fVarsForOpt = NULL;
632 //printf("Changing the variables for cut optimization\n");
633 }
634
635 if(nVars==0){//!=fnVars) {
636 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
637 return;
638 }
639
640 fnVarsForOpt = 0;
641 fVarsForOpt = new Bool_t[fnVars];
642 for(Int_t iv=0; iv<fnVars; iv++) {
643 fVarsForOpt[iv]=forOpt[iv];
644 if(fVarsForOpt[iv]) fnVarsForOpt++;
645 }
646
647 return;
648}
649
650//---------------------------------------------------------------------------
651void AliRDHFCuts::SetUseCentrality(Int_t flag) {
652 //
653 // set centrality estimator
654 //
655 fUseCentrality=flag;
656 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
657
658 return;
659}
660
661
662//---------------------------------------------------------------------------
663void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
664 //
665 // store the cuts
666 //
667 if(nVars!=fnVars) {
668 printf("Wrong number of variables: it has to be %d\n",fnVars);
669 AliFatal("exiting");
670 }
671 if(nPtBins!=fnPtBins) {
672 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
673 AliFatal("exiting");
674 }
675
676 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
677
678
679 for(Int_t iv=0; iv<fnVars; iv++) {
680
681 for(Int_t ib=0; ib<fnPtBins; ib++) {
682
683 //check
684 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
685 cout<<"Overflow, exit..."<<endl;
686 return;
687 }
688
689 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
690
691 }
692 }
693 return;
694}
695//---------------------------------------------------------------------------
696void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
697 //
698 // store the cuts
699 //
700 if(glIndex != fGlobalIndex){
701 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
702 AliFatal("exiting");
703 }
704 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
705
706 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
707 fCutsRD[iGl] = cutsRDGlob[iGl];
708 }
709 return;
710}
711//---------------------------------------------------------------------------
712void AliRDHFCuts::PrintAll() const {
713 //
714 // print all cuts values
715 //
716
717 printf("Minimum vtx type %d\n",fMinVtxType);
718 printf("Minimum vtx contr %d\n",fMinVtxContr);
719 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
720 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
721 printf("Use PID %d\n",(Int_t)fUsePID);
722 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
4bdc0315 723 printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
0c22e2ac 724 printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
725 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
726 if(fOptPileup==1) printf(" -- Reject pileup event");
727 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
728 if(fUseCentrality>0) {
729 TString estimator="";
730 if(fUseCentrality==1) estimator = "V0";
731 if(fUseCentrality==2) estimator = "Tracks";
732 if(fUseCentrality==3) estimator = "Tracklets";
733 if(fUseCentrality==4) estimator = "SPD clusters outer";
734 printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
735 }
736 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
737
738 if(fVarNames){
739 cout<<"Array of variables"<<endl;
740 for(Int_t iv=0;iv<fnVars;iv++){
741 cout<<fVarNames[iv]<<"\t";
742 }
743 cout<<endl;
744 }
745 if(fVarsForOpt){
746 cout<<"Array of optimization"<<endl;
747 for(Int_t iv=0;iv<fnVars;iv++){
748 cout<<fVarsForOpt[iv]<<"\t";
749 }
750 cout<<endl;
751 }
752 if(fIsUpperCut){
753 cout<<"Array of upper/lower cut"<<endl;
754 for(Int_t iv=0;iv<fnVars;iv++){
755 cout<<fIsUpperCut[iv]<<"\t";
756 }
757 cout<<endl;
758 }
759 if(fPtBinLimits){
760 cout<<"Array of ptbin limits"<<endl;
761 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
762 cout<<fPtBinLimits[ib]<<"\t";
763 }
764 cout<<endl;
765 }
766 if(fCutsRD){
767 cout<<"Matrix of cuts"<<endl;
768 for(Int_t iv=0;iv<fnVars;iv++){
769 for(Int_t ib=0;ib<fnPtBins;ib++){
770 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
771 }
772 cout<<endl;
773 }
774 cout<<endl;
775 }
776 return;
777}
778
779//--------------------------------------------------------------------------
780void AliRDHFCuts::PrintTrigger() const{
781
782 cout<<" Trigger selection pattern: ";
783 if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
784 if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
785 if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
786 if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
787 if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
788 if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
789 if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
790 if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
791 if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
792 if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
793 cout << endl<< endl;
794
795}
796
797//---------------------------------------------------------------------------
798void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
799 //
800 // get the cuts
801 //
802
803 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
804
805
806 Int_t iv,ib;
807 if(!cutsRD) {
808 //cout<<"Initialization..."<<endl;
809 cutsRD=new Float_t*[fnVars];
810 for(iv=0; iv<fnVars; iv++) {
811 cutsRD[iv] = new Float_t[fnPtBins];
812 }
813 }
814
815 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
816 GetVarPtIndex(iGlobal,iv,ib);
817 cutsRD[iv][ib] = fCutsRD[iGlobal];
818 }
819
820 return;
821}
822
823//---------------------------------------------------------------------------
824Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
825 //
826 // give the global index from variable and pt bin
827 //
828 return iPtBin*fnVars+iVar;
829}
830
831//---------------------------------------------------------------------------
832void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
833 //
834 //give the index of the variable and of the pt bin from the global index
835 //
836 iPtBin=(Int_t)iGlob/fnVars;
837 iVar=iGlob%fnVars;
838
839 return;
840}
841
842//---------------------------------------------------------------------------
843Int_t AliRDHFCuts::PtBin(Double_t pt) const {
844 //
845 //give the pt bin where the pt lies.
846 //
847 Int_t ptbin=-1;
848 if(pt<fPtBinLimits[0])return ptbin;
849 for (Int_t i=0;i<fnPtBins;i++){
850 if(pt<fPtBinLimits[i+1]) {
851 ptbin=i;
852 break;
853 }
854 }
855 return ptbin;
856}
857//-------------------------------------------------------------------
858Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
859 //
860 // Give the value of cut set for the variable iVar and the pt bin iPtBin
861 //
862 if(!fCutsRD){
863 cout<<"Cuts not iniziaisez yet"<<endl;
864 return 0;
865 }
866 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
867}
868//-------------------------------------------------------------------
869Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
870 //
871 // Get centrality percentile
872 //
873
874 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
875 if(mcArray) {fUseAOD049=kFALSE;}
876
877 AliAODHeader *header=aodEvent->GetHeader();
878 AliCentrality *centrality=header->GetCentralityP();
879 Float_t cent=-999.;
880 Bool_t isSelRun=kFALSE;
881 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
882 if(!centrality) return cent;
883 else{
884 if (estimator==kCentV0M){
885 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
886 if(cent<0){
887 Int_t quality = centrality->GetQuality();
888 if(quality<=1){
889 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
890 }else{
891 Int_t runnum=aodEvent->GetRunNumber();
892 for(Int_t ir=0;ir<5;ir++){
893 if(runnum==selRun[ir]){
894 isSelRun=kTRUE;
895 break;
896 }
897 }
898 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
899 }
900 }
901
902 //temporary fix for AOD049 outliers
903 if(fUseAOD049&&cent>=0){
904 Float_t v0=0;
905 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
906 v0+=aodV0->GetMTotV0A();
907 v0+=aodV0->GetMTotV0C();
908 if(cent==0&&v0<19500)return -1;//filtering issue
909 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
910 Float_t val= 1.30552 + 0.147931 * v0;
911 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};
912 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
913 }
914 }
915 else {
916 if (estimator==kCentTRK) {
917 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
918 if(cent<0){
919 Int_t quality = centrality->GetQuality();
920 if(quality<=1){
921 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
922 }else{
923 Int_t runnum=aodEvent->GetRunNumber();
924 for(Int_t ir=0;ir<5;ir++){
925 if(runnum==selRun[ir]){
926 isSelRun=kTRUE;
927 break;
928 }
929 }
930 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
931 }
932 }
933 }
934 else{
935 if (estimator==kCentTKL){
936 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
937 if(cent<0){
938 Int_t quality = centrality->GetQuality();
939 if(quality<=1){
940 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
941 }else{
942 Int_t runnum=aodEvent->GetRunNumber();
943 for(Int_t ir=0;ir<5;ir++){
944 if(runnum==selRun[ir]){
945 isSelRun=kTRUE;
946 break;
947 }
948 }
949 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
950 }
951 }
952 }
953 else{
954 if (estimator==kCentCL1){
955 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
956 if(cent<0){
957 Int_t quality = centrality->GetQuality();
958 if(quality<=1){
959 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
960 }else{
961 Int_t runnum=aodEvent->GetRunNumber();
962 for(Int_t ir=0;ir<5;ir++){
963 if(runnum==selRun[ir]){
964 isSelRun=kTRUE;
965 break;
966 }
967 }
968 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
969 }
970 }
971 }
972 else {
973 AliWarning("Centrality estimator not valid");
974
975 }
976 }
977 }
978 }
979 }
980 return cent;
981}
982//-------------------------------------------------------------------
983Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
984 //
985 // Compare two cuts objects
986 //
987
988 Bool_t areEqual=kTRUE;
989
990 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
991
992 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
993
994 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
995
996 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
997
998 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
999
1000 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1001 if(fTrackCuts){
1002 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1003
1004 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1005
1006 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1007
1008 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;}
1009 }
1010
1011 if(fCutsRD) {
1012 for(Int_t iv=0;iv<fnVars;iv++) {
1013 for(Int_t ib=0;ib<fnPtBins;ib++) {
1014 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1015 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1016 areEqual=kFALSE;
1017 }
1018 }
1019 }
1020 }
1021
1022 return areEqual;
1023}
1024//---------------------------------------------------------------------------
1025void AliRDHFCuts::MakeTable() const {
1026 //
1027 // print cuts values in table format
1028 //
1029
1030 TString ptString = "pT range";
1031 if(fVarNames && fPtBinLimits && fCutsRD){
1032 TString firstLine(Form("* %-15s",ptString.Data()));
1033 for (Int_t ivar=0; ivar<fnVars; ivar++){
1034 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1035 if (ivar == fnVars){
1036 firstLine+="*\n";
1037 }
1038 }
1039 Printf("%s",firstLine.Data());
1040
1041 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1042 TString line;
1043 if (ipt==fnPtBins-1){
1044 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1045 }
1046 else{
1047 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1048 }
1049 for (Int_t ivar=0; ivar<fnVars; ivar++){
1050 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1051 }
1052 Printf("%s",line.Data());
1053 }
1054
1055 }
1056
1057
1058 return;
1059}
1060//--------------------------------------------------------------------------
1061Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1062 AliAODEvent *aod) const
1063{
1064 //
1065 // Recalculate primary vertex without daughters
1066 //
1067
1068 if(!aod) {
1069 AliError("Can not remove daughters from vertex without AOD event");
1070 return 0;
1071 }
1072
1073 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1074 if(!recvtx){
1075 AliDebug(2,"Removal of daughter tracks failed");
1076 return kFALSE;
1077 }
1078
1079
1080 //set recalculed primary vertex
1081 d->SetOwnPrimaryVtx(recvtx);
1082 delete recvtx;
1083
1084 return kTRUE;
1085}
1086//--------------------------------------------------------------------------
1087Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1088{
1089 //
1090 // Recalculate primary vertex without daughters
1091 //
1092
1093 if(!aod) {
1094 AliError("Can not get MC vertex without AOD event");
1095 return kFALSE;
1096 }
1097
1098 // load MC header
1099 AliAODMCHeader *mcHeader =
1100 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1101 if(!mcHeader) {
1102 AliError("Can not get MC vertex without AODMCHeader event");
1103 return kFALSE;
1104 }
1105 Double_t pos[3];
1106 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1107 mcHeader->GetVertex(pos);
1108 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1109
1110 if(!recvtx){
1111 AliDebug(2,"Removal of daughter tracks failed");
1112 return kFALSE;
1113 }
1114
1115 //set recalculed primary vertex
1116 d->SetOwnPrimaryVtx(recvtx);
1117
1118 d->RecalculateImpPars(recvtx,aod);
1119
1120 delete recvtx;
1121
1122 return kTRUE;
1123}
1124//--------------------------------------------------------------------------
1125void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1126 AliAODEvent *aod,
1127 AliAODVertex *origownvtx) const
1128{
1129 //
1130 // Clean-up own primary vertex if needed
1131 //
1132
1133 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1134 d->UnsetOwnPrimaryVtx();
1135 if(origownvtx) {
1136 d->SetOwnPrimaryVtx(origownvtx);
1137 delete origownvtx; origownvtx=NULL;
1138 }
1139 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1140 } else {
1141 if(origownvtx) {
1142 delete origownvtx; origownvtx=NULL;
1143 }
1144 }
1145 return;
1146}
1147//--------------------------------------------------------------------------
1148Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1149{
1150 //
1151 // Checks if this candidate is matched to MC signal
1152 //
1153
1154 if(!aod) return kFALSE;
1155
1156 // get the MC array
1157 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1158
1159 if(!mcArray) return kFALSE;
1160
1161 // try to match
1162 Int_t label = d->MatchToMC(pdg,mcArray);
1163
1164 if(label>=0) {
1165 //printf("MATCH!\n");
1166 return kTRUE;
1167 }
1168
1169 return kFALSE;
1170}
1171
1172
4bdc0315 1173//--------------------------------------------------------------------------
1174Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1175 // recompute event primary vertex from AOD tracks
1176
1177 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1178 vertexer->SetITSMode();
1179 vertexer->SetMinClusters(3);
1180
1181 AliAODVertex* pvtx=event->GetPrimaryVertex();
1182 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1183 Float_t diamondcovxy[3];
1184 event->GetDiamondCovXY(diamondcovxy);
1185 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1186 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1187 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1188 vertexer->SetVtxStart(diamond);
1189 delete diamond; diamond=NULL;
1190 }
1191
1192 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1193 if(!vertexESD) return kFALSE;
1194 if(vertexESD->GetNContributors()<=0) {
1195 //AliDebug(2,"vertexing failed");
1196 delete vertexESD; vertexESD=NULL;
1197 return kFALSE;
1198 }
1199 delete vertexer; vertexer=NULL;
1200
1201 // convert to AliAODVertex
1202 Double_t pos[3],cov[6],chi2perNDF;
1203 vertexESD->GetXYZ(pos); // position
1204 vertexESD->GetCovMatrix(cov); //covariance matrix
1205 chi2perNDF = vertexESD->GetChi2toNDF();
1206 delete vertexESD; vertexESD=NULL;
1207
1208 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1209 pvtx->SetChi2perNDF(chi2perNDF);
1210 pvtx->SetCovMatrix(cov);
1211
1212 return kTRUE;
1213}