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