1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // Base class for cuts on AOD reconstructed heavy-flavour decay
22 // Author: A.Dainese, andrea.dainese@pd.infn.it
23 /////////////////////////////////////////////////////////////
24 #include <Riostream.h>
26 #include "AliVEvent.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliVVertex.h"
30 #include "AliESDVertex.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"
41 #include "AliVertexerTracks.h"
42 #include "AliRDHFCuts.h"
43 #include "AliAnalysisManager.h"
44 #include "AliInputEventHandler.h"
45 #include "AliPIDResponse.h"
53 //--------------------------------------------------------------------------
54 AliRDHFCuts::AliRDHFCuts(const Char_t* name, const Char_t* title) :
55 AliAnalysisCuts(name,title),
60 fMinSPDMultiplicity(0),
61 fTriggerMask(AliVEvent::kAnyINT),
62 fUseOnlyOneTrigger(kFALSE),
79 fRemoveDaughtersFromPrimary(kFALSE),
80 fRecomputePrimVertex(kFALSE),
82 fUsePhysicsSelection(kTRUE),
94 fKeepSignalMC(kFALSE),
95 fIsCandTrackSPDFirst(kFALSE),
96 fMaxPtCandTrackSPDFirst(0.),
97 fApplySPDDeadPbPb2011(kFALSE),
98 fRemoveTrackletOutliers(kFALSE),
101 fUseTrackSelectionWithFilterBits(kTRUE),
102 fUseCentrFlatteningInMC(kFALSE),
106 // Default Constructor
108 fTriggerClass[0]="CINT1"; fTriggerClass[1]="";
110 //--------------------------------------------------------------------------
111 AliRDHFCuts::AliRDHFCuts(const AliRDHFCuts &source) :
112 AliAnalysisCuts(source),
113 fMinVtxType(source.fMinVtxType),
114 fMinVtxContr(source.fMinVtxContr),
115 fMaxVtxRedChi2(source.fMaxVtxRedChi2),
116 fMaxVtxZ(source.fMaxVtxZ),
117 fMinSPDMultiplicity(source.fMinSPDMultiplicity),
118 fTriggerMask(source.fTriggerMask),
119 fUseOnlyOneTrigger(source.fUseOnlyOneTrigger),
122 fnPtBins(source.fnPtBins),
123 fnPtBinLimits(source.fnPtBinLimits),
125 fnVars(source.fnVars),
127 fnVarsForOpt(source.fnVarsForOpt),
129 fGlobalIndex(source.fGlobalIndex),
132 fUsePID(source.fUsePID),
133 fUseAOD049(source.fUseAOD049),
135 fWhyRejection(source.fWhyRejection),
136 fEvRejectionBits(source.fEvRejectionBits),
137 fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
138 fRecomputePrimVertex(source.fRecomputePrimVertex),
139 fUseMCVertex(source.fUseMCVertex),
140 fUsePhysicsSelection(source.fUsePhysicsSelection),
141 fOptPileup(source.fOptPileup),
142 fMinContrPileup(source.fMinContrPileup),
143 fMinDzPileup(source.fMinDzPileup),
144 fUseCentrality(source.fUseCentrality),
145 fMinCentrality(source.fMinCentrality),
146 fMaxCentrality(source.fMaxCentrality),
147 fFixRefs(source.fFixRefs),
148 fIsSelectedCuts(source.fIsSelectedCuts),
149 fIsSelectedPID(source.fIsSelectedPID),
150 fMinPtCand(source.fMinPtCand),
151 fMaxPtCand(source.fMaxPtCand),
152 fKeepSignalMC(source.fKeepSignalMC),
153 fIsCandTrackSPDFirst(source.fIsCandTrackSPDFirst),
154 fMaxPtCandTrackSPDFirst(source.fMaxPtCandTrackSPDFirst),
155 fApplySPDDeadPbPb2011(source.fApplySPDDeadPbPb2011),
156 fRemoveTrackletOutliers(source.fRemoveTrackletOutliers),
157 fCutOnzVertexSPD(source.fCutOnzVertexSPD),
158 fKinkReject(source.fKinkReject),
159 fUseTrackSelectionWithFilterBits(source.fUseTrackSelectionWithFilterBits),
160 fUseCentrFlatteningInMC(source.fUseCentrFlatteningInMC),
166 cout<<"Copy constructor"<<endl;
167 fTriggerClass[0] = source.fTriggerClass[0];
168 fTriggerClass[1] = source.fTriggerClass[1];
169 if(source.GetTrackCuts()) AddTrackCuts(source.GetTrackCuts());
170 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
171 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
172 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
173 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
174 if(source.fPidHF) SetPidHF(source.fPidHF);
175 if(source.fHistCentrDistr) fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
179 //--------------------------------------------------------------------------
180 AliRDHFCuts &AliRDHFCuts::operator=(const AliRDHFCuts &source)
183 // assignment operator
185 if(&source == this) return *this;
187 AliAnalysisCuts::operator=(source);
189 fMinVtxType=source.fMinVtxType;
190 fMinVtxContr=source.fMinVtxContr;
191 fMaxVtxRedChi2=source.fMaxVtxRedChi2;
192 fMaxVtxZ=source.fMaxVtxZ;
193 fMinSPDMultiplicity=source.fMinSPDMultiplicity;
194 fTriggerMask=source.fTriggerMask;
195 fUseOnlyOneTrigger=source.fUseOnlyOneTrigger;
196 fTriggerClass[0]=source.fTriggerClass[0];
197 fTriggerClass[1]=source.fTriggerClass[1];
198 fnPtBins=source.fnPtBins;
199 fnPtBinLimits=source.fnPtBinLimits;
200 fnVars=source.fnVars;
201 fGlobalIndex=source.fGlobalIndex;
202 fnVarsForOpt=source.fnVarsForOpt;
203 fUsePID=source.fUsePID;
204 fUseAOD049=source.fUseAOD049;
205 if(fPidHF) delete fPidHF;
206 fPidHF=new AliAODPidHF(*(source.GetPidHF()));
207 fWhyRejection=source.fWhyRejection;
208 fEvRejectionBits=source.fEvRejectionBits;
209 fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
210 fRecomputePrimVertex=source.fRecomputePrimVertex;
211 fUseMCVertex=source.fUseMCVertex;
212 fUsePhysicsSelection=source.fUsePhysicsSelection;
213 fOptPileup=source.fOptPileup;
214 fMinContrPileup=source.fMinContrPileup;
215 fMinDzPileup=source.fMinDzPileup;
216 fUseCentrality=source.fUseCentrality;
217 fMinCentrality=source.fMinCentrality;
218 fMaxCentrality=source.fMaxCentrality;
219 fFixRefs=source.fFixRefs;
220 fIsSelectedCuts=source.fIsSelectedCuts;
221 fIsSelectedPID=source.fIsSelectedPID;
222 fMinPtCand=source.fMinPtCand;
223 fMaxPtCand=source.fMaxPtCand;
224 fKeepSignalMC=source.fKeepSignalMC;
225 fIsCandTrackSPDFirst=source.fIsCandTrackSPDFirst;
226 fMaxPtCandTrackSPDFirst=source.fMaxPtCandTrackSPDFirst;
227 fApplySPDDeadPbPb2011=source.fApplySPDDeadPbPb2011;
228 fRemoveTrackletOutliers=source.fRemoveTrackletOutliers;
229 fCutOnzVertexSPD=source.fCutOnzVertexSPD;
230 fKinkReject=source.fKinkReject;
231 fUseTrackSelectionWithFilterBits=source.fUseTrackSelectionWithFilterBits;
232 if(fHistCentrDistr) delete fHistCentrDistr;
233 fUseCentrFlatteningInMC=source.fUseCentrFlatteningInMC;
234 if(source.fHistCentrDistr)fHistCentrDistr=(TH1F*)(source.fHistCentrDistr->Clone());
236 if(source.GetTrackCuts()) {delete fTrackCuts; fTrackCuts=new AliESDtrackCuts(*(source.GetTrackCuts()));}
237 if(source.fPtBinLimits) SetPtBins(source.fnPtBinLimits,source.fPtBinLimits);
238 if(source.fVarNames) SetVarNames(source.fnVars,source.fVarNames,source.fIsUpperCut);
239 if(source.fCutsRD) SetCuts(source.fGlobalIndex,source.fCutsRD);
240 if(source.fVarsForOpt) SetVarsForOpt(source.fnVarsForOpt,source.fVarsForOpt);
245 //--------------------------------------------------------------------------
246 AliRDHFCuts::~AliRDHFCuts() {
248 // Default Destructor
250 if(fPtBinLimits) {delete [] fPtBinLimits; fPtBinLimits=0;}
251 if(fVarNames) {delete [] fVarNames; fVarNames=0;}
252 if(fVarsForOpt) {delete [] fVarsForOpt; fVarsForOpt=0;}
257 if(fIsUpperCut) {delete [] fIsUpperCut; fIsUpperCut=0;}
262 if(fHistCentrDistr)delete fHistCentrDistr;
264 //---------------------------------------------------------------------------
265 Int_t AliRDHFCuts::IsEventSelectedInCentrality(AliVEvent *event) {
267 // Centrality selection
269 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid){
270 AliWarning("Centrality estimator not valid");
273 Float_t centvalue=GetCentrality((AliAODEvent*)event);
274 if (centvalue<-998.){//-999 if no centralityP
277 if (centvalue<fMinCentrality || centvalue>fMaxCentrality){
280 // selection to flatten the centrality distribution
282 if(!IsEventSelectedForCentrFlattening(centvalue))return 4;
290 //-------------------------------------------------
291 void AliRDHFCuts::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
292 // set the histo for centrality flattening
293 // the centrality is flatten in the range minCentr,maxCentr
294 // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
295 // positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
296 // negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
297 // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
299 if(maxCentr<minCentr){
300 AliWarning("AliRDHFCuts::Wrong centralities values while setting the histogram for centrality flattening");
303 if(fHistCentrDistr)delete fHistCentrDistr;
304 fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
305 fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
306 Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
307 Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
308 fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
309 Double_t ref=0.,bincont=0.,binrefwidth=1.;
311 if(TMath::Abs(centrRef)<0.0001){
312 binref=fHistCentrDistr->GetMinimumBin();
313 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
314 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
316 else if(centrRef>0.){
317 binref=h->FindBin(centrRef);
318 if(binref<1||binref>h->GetNbinsX()){
319 AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
321 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
322 ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
325 if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
326 binref=fHistCentrDistr->GetMaximumBin();
327 binrefwidth=fHistCentrDistr->GetBinWidth(binref);
328 ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
331 for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
332 if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
333 bincont=h->GetBinContent(j);
334 fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
335 fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
338 h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
342 fHistCentrDistr->SetBinContent(0,switchTRand);
347 //-------------------------------------------------
348 Bool_t AliRDHFCuts::IsEventSelectedForCentrFlattening(Float_t centvalue){
350 // Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
351 // Can be faster if it was required that fHistCentrDistr covers
352 // exactly the desired centrality range (e.g. part of the lines below should be done during the
353 // setting of the histo) and TH1::SetMinimum called
356 if(!fHistCentrDistr) return kTRUE;
357 // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
358 // if(maxbin>fHistCentrDistr->GetNbinsX()){
359 // AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
362 Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
363 Double_t bincont=fHistCentrDistr->GetBinContent(bin);
364 Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
366 if(fHistCentrDistr->GetBinContent(0)<-0.9999){
367 if(gRandom->Uniform(1.)<bincont)return kTRUE;
371 if(centDigits*100.<bincont)return kTRUE;
375 //---------------------------------------------------------------------------
376 void AliRDHFCuts::SetupPID(AliVEvent *event) {
377 // Set the PID response object in the AliAODPidHF
378 // in case of old PID sets the TPC dE/dx BB parameterization
381 if(fPidHF->GetPidResponse()==0x0){
382 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
383 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
384 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
385 fPidHF->SetPidResponse(pidResp);
387 if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
388 if(fPidHF->GetOldPid()) {
391 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
392 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
394 // pp, from LHC10d onwards
395 if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
396 event->GetRunNumber()>139517) fPidHF->SetOnePad(kTRUE);
397 // pp, 2011 low energy run
398 if((event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860)){
399 fPidHF->SetppLowEn2011(kTRUE);
400 fPidHF->SetOnePad(kFALSE);
403 if(event->GetRunNumber()>=136851 && event->GetRunNumber()<=139517) fPidHF->SetPbPb(kTRUE);
405 if(isMC) fPidHF->SetMC(kTRUE);
406 if(isMC && (event->GetRunNumber()>=146686 && event->GetRunNumber()<=146860))
407 fPidHF->SetMClowenpp2011(kTRUE);
408 fPidHF->SetBetheBloch();
410 // check that AliPIDResponse object was properly set in case of using OADB
411 if(fPidHF->GetPidResponse()==0x0) AliFatal("AliPIDResponse object not set");
415 //---------------------------------------------------------------------------
416 Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
420 //if(fTriggerMask && event->GetTriggerMask()!=fTriggerMask) return kFALSE;
428 if(fRecomputePrimVertex){
429 Bool_t vertOK= RecomputePrimaryVertex((AliAODEvent*)event);
438 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)event)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
439 if(mcArray) {isMC=kTRUE;fUseAOD049=kFALSE;}
445 TString firedTriggerClasses=((AliAODEvent*)event)->GetFiredTriggerClasses();
446 // don't do for MC and for PbPb 2010 data
447 if(!isMC && (event->GetRunNumber()<136851 || event->GetRunNumber()>139517)) {
448 if(!firedTriggerClasses.Contains(fTriggerClass[0].Data()) &&
449 (fTriggerClass[1].CompareTo("")==0 || !firedTriggerClasses.Contains(fTriggerClass[1].Data())) ) {
451 fEvRejectionBits+=1<<kNotSelTrigger;
456 // TEMPORARY FIX FOR GetEvent
457 Int_t nTracks=((AliAODEvent*)event)->GetNTracks();
458 for(Int_t itr=0; itr<nTracks; itr++){
459 AliAODTrack* tr=(AliAODTrack*)((AliAODEvent*)event)->GetTrack(itr);
460 tr->SetAODEvent((AliAODEvent*)event);
463 // TEMPORARY FIX FOR REFERENCES
464 // Fix references to daughter tracks
466 // AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
467 // fixer->FixReferences((AliAODEvent*)event);
473 // physics selection requirements
474 if(fUsePhysicsSelection){
475 Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerMask);
477 if(accept) fWhyRejection=7;
478 fEvRejectionBits+=1<<kPhysicsSelection;
481 if(fUseOnlyOneTrigger){
482 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()!=fTriggerMask){
483 if(accept) fWhyRejection=7;
484 fEvRejectionBits+=1<<kPhysicsSelection;
491 // vertex requirements
493 const AliVVertex *vertex = event->GetPrimaryVertex();
497 fEvRejectionBits+=1<<kNoVertex;
499 TString title=vertex->GetTitle();
500 if(title.Contains("Z") && fMinVtxType>1){
502 fEvRejectionBits+=1<<kNoVertex;
504 else if(title.Contains("3D") && fMinVtxType>2){
506 fEvRejectionBits+=1<<kNoVertex;
508 if(vertex->GetNContributors()<fMinVtxContr){
510 fEvRejectionBits+=1<<kTooFewVtxContrib;
512 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
513 fEvRejectionBits+=1<<kZVtxOutFid;
514 if(accept) fWhyRejection=6;
519 if(fCutOnzVertexSPD>0){
520 const AliVVertex *vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
521 if(!vSPD || (vSPD && vSPD->GetNContributors()<fMinVtxContr)){
523 fEvRejectionBits+=1<<kBadSPDVertex;
525 if(fCutOnzVertexSPD==1 && TMath::Abs(vSPD->GetZ())>12.) {
526 fEvRejectionBits+=1<<kZVtxSPDOutFid;
527 if(accept) fWhyRejection=6;
530 if(fCutOnzVertexSPD==2 && vertex){
531 if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) {
532 fEvRejectionBits+=1<<kZVtxSPDOutFid;
533 if(accept) fWhyRejection=6;
541 if(fOptPileup==kRejectPileupEvent){
542 Int_t cutc=(Int_t)fMinContrPileup;
543 Double_t cutz=(Double_t)fMinDzPileup;
544 if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) {
545 if(accept) fWhyRejection=1;
546 fEvRejectionBits+=1<<kPileupSPD;
551 // centrality selection
552 if (fUseCentrality!=kCentOff) {
553 Int_t rejection=IsEventSelectedInCentrality(event);
554 Bool_t okCent=kFALSE;
555 if(rejection==0) okCent=kTRUE;
556 if(isMC && rejection==4 && !fUseCentrFlatteningInMC) okCent=kTRUE;
558 if(accept) fWhyRejection=rejection;
559 if(fWhyRejection==4)fEvRejectionBits+=1<<kCentralityFlattening;
560 else fEvRejectionBits+=1<<kOutsideCentrality;
566 // PbPb2011 outliers in tracklets vs. VZERO
567 if(event->GetRunNumber()>=167693 && event->GetRunNumber()<=170593){
568 if(fRemoveTrackletOutliers){
569 Double_t v0cent=GetCentrality((AliAODEvent*)event,kCentV0M);
570 Double_t ntracklets=((AliAODEvent*)event)->GetTracklets()->GetNumberOfTracklets();
571 Double_t cutval=60.-0.08*ntracklets+1./50000.*ntracklets*ntracklets;
572 if(ntracklets<1000. && v0cent<cutval){
573 if(accept) fWhyRejection=2;
574 fEvRejectionBits+=1<<kOutsideCentrality;
582 //---------------------------------------------------------------------------
583 Bool_t AliRDHFCuts::AreDaughtersSelected(AliAODRecoDecayHF *d) const {
585 // Daughter tracks selection
587 if(!fTrackCuts) return kTRUE;
589 Int_t ndaughters = d->GetNDaughters();
590 AliAODVertex *vAOD = d->GetPrimaryVtx();
591 Double_t pos[3],cov[6];
593 vAOD->GetCovarianceMatrix(cov);
594 const AliESDVertex vESD(pos,cov,100.,100);
598 for(Int_t idg=0; idg<ndaughters; idg++) {
599 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(idg);
600 if(!dgTrack) {retval = kFALSE; continue;}
601 //printf("charge %d\n",dgTrack->Charge());
602 if(dgTrack->Charge()==0) continue; // it's not a track, but a V0
604 if(fIsCandTrackSPDFirst && d->Pt()<fMaxPtCandTrackSPDFirst)
605 { if(!dgTrack->HasPointOnITSLayer(0)) { retval = kFALSE; continue; } }
607 if(!IsDaughterSelected(dgTrack,&vESD,fTrackCuts)) retval = kFALSE;
612 //---------------------------------------------------------------------------
613 Bool_t AliRDHFCuts::IsDaughterSelected(AliAODTrack *track,const AliESDVertex *primary,AliESDtrackCuts *cuts) const {
615 // Convert to ESDtrack, relate to vertex and check cuts
617 if(!cuts) return kTRUE;
619 if(cuts->GetFlagCutTOFdistance()) cuts->SetFlagCutTOFdistance(kFALSE);
623 // convert to ESD track here
624 AliESDtrack esdTrack(track);
625 // set the TPC cluster info
626 esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
627 esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
628 esdTrack.SetTPCPointsF(track->GetTPCNclsF());
629 // needed to calculate the impact parameters
630 esdTrack.RelateToVertex(primary,0.,3.);
632 if(!cuts->IsSelected(&esdTrack)) retval = kFALSE;
636 AliAODVertex *maybeKink=track->GetProdVertex();
637 if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
640 if(fOptPileup==kRejectTracksFromPileupVertex){
642 // we need either to have here the AOD Event,
643 // or to have the pileup vertex object
646 if(fApplySPDDeadPbPb2011){
647 Bool_t deadSPDLay1PbPb2011[20][4]={
648 {kTRUE,kTRUE,kTRUE,kTRUE},
649 {kTRUE,kTRUE,kTRUE,kTRUE},
650 {kTRUE,kTRUE,kTRUE,kTRUE},
651 {kTRUE,kTRUE,kTRUE,kTRUE},
652 {kTRUE,kTRUE,kTRUE,kTRUE},
653 {kFALSE,kFALSE,kTRUE,kTRUE},
654 {kTRUE,kTRUE,kFALSE,kFALSE},
655 {kTRUE,kTRUE,kTRUE,kTRUE},
656 {kFALSE,kFALSE,kTRUE,kTRUE},
657 {kTRUE,kTRUE,kTRUE,kTRUE},
658 {kTRUE,kTRUE,kFALSE,kFALSE},
659 {kTRUE,kTRUE,kTRUE,kTRUE},
660 {kFALSE,kFALSE,kFALSE,kFALSE},
661 {kFALSE,kFALSE,kTRUE,kTRUE},
662 {kFALSE,kFALSE,kFALSE,kFALSE},
663 {kFALSE,kFALSE,kFALSE,kFALSE},
664 {kTRUE,kTRUE,kTRUE,kTRUE},
665 {kTRUE,kTRUE,kFALSE,kFALSE},
666 {kFALSE,kFALSE,kFALSE,kFALSE},
667 {kFALSE,kFALSE,kFALSE,kFALSE}
669 Bool_t deadSPDLay2PbPb2011[40][4]={
670 {kTRUE,kTRUE,kTRUE,kTRUE},
671 {kTRUE,kTRUE,kTRUE,kTRUE},
672 {kTRUE,kTRUE,kTRUE,kTRUE},
673 {kTRUE,kTRUE,kTRUE,kTRUE},
674 {kTRUE,kTRUE,kTRUE,kTRUE},
675 {kTRUE,kTRUE,kTRUE,kTRUE},
676 {kTRUE,kTRUE,kTRUE,kTRUE},
677 {kTRUE,kTRUE,kTRUE,kTRUE},
678 {kTRUE,kTRUE,kTRUE,kTRUE},
679 {kTRUE,kTRUE,kTRUE,kTRUE},
680 {kTRUE,kTRUE,kTRUE,kTRUE},
681 {kTRUE,kTRUE,kTRUE,kTRUE},
682 {kFALSE,kFALSE,kFALSE,kFALSE},
683 {kFALSE,kFALSE,kTRUE,kTRUE},
684 {kTRUE,kTRUE,kTRUE,kTRUE},
685 {kTRUE,kTRUE,kTRUE,kTRUE},
686 {kTRUE,kTRUE,kFALSE,kFALSE},
687 {kTRUE,kTRUE,kTRUE,kTRUE},
688 {kTRUE,kTRUE,kTRUE,kTRUE},
689 {kTRUE,kTRUE,kTRUE,kTRUE},
690 {kFALSE,kFALSE,kFALSE,kFALSE},
691 {kFALSE,kFALSE,kFALSE,kFALSE},
692 {kTRUE,kTRUE,kTRUE,kTRUE},
693 {kTRUE,kTRUE,kTRUE,kTRUE},
694 {kFALSE,kFALSE,kFALSE,kFALSE},
695 {kFALSE,kFALSE,kFALSE,kFALSE},
696 {kTRUE,kTRUE,kTRUE,kTRUE},
697 {kTRUE,kTRUE,kTRUE,kTRUE},
698 {kFALSE,kFALSE,kFALSE,kFALSE},
699 {kFALSE,kFALSE,kFALSE,kFALSE},
700 {kFALSE,kFALSE,kFALSE,kFALSE},
701 {kFALSE,kFALSE,kFALSE,kFALSE},
702 {kTRUE,kTRUE,kTRUE,kTRUE},
703 {kTRUE,kTRUE,kTRUE,kTRUE},
704 {kTRUE,kTRUE,kFALSE,kFALSE},
705 {kTRUE,kTRUE,kTRUE,kTRUE},
706 {kFALSE,kFALSE,kFALSE,kFALSE},
707 {kFALSE,kFALSE,kFALSE,kFALSE},
708 {kFALSE,kFALSE,kFALSE,kFALSE},
709 {kFALSE,kFALSE,kFALSE,kFALSE}
711 Double_t xyz1[3],xyz2[3];
712 esdTrack.GetXYZAt(3.9,0.,xyz1);
713 esdTrack.GetXYZAt(7.6,0.,xyz2);
714 Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
715 if(phi1<0) phi1+=2*TMath::Pi();
716 Int_t lad1=(Int_t)(phi1/(2.*TMath::Pi()/20.));
717 Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
718 if(phi2<0) phi2+=2*TMath::Pi();
719 Int_t lad2=(Int_t)(phi2/(2.*TMath::Pi()/40.));
720 Int_t mod1=(Int_t)((xyz1[2]+14)/7.);
721 Int_t mod2=(Int_t)((xyz2[2]+14)/7.);
722 Bool_t lay1ok=kFALSE;
723 if(mod1>=0 && mod1<4 && lad1<20){
724 lay1ok=deadSPDLay1PbPb2011[lad1][mod1];
726 Bool_t lay2ok=kFALSE;
727 if(mod2>=0 && mod2<4 && lad2<40){
728 lay2ok=deadSPDLay2PbPb2011[lad2][mod2];
730 if(!lay1ok && !lay2ok) retval=kFALSE;
735 //---------------------------------------------------------------------------
736 void AliRDHFCuts::SetPtBins(Int_t nPtBinLimits,Float_t *ptBinLimits) {
740 delete [] fPtBinLimits;
742 printf("Changing the pt bins\n");
745 if(nPtBinLimits != fnPtBins+1){
746 cout<<"Warning: ptBinLimits dimention "<<nPtBinLimits<<" != nPtBins+1 ("<<fnPtBins+1<<")\nSetting nPtBins to "<<nPtBinLimits-1<<endl;
747 SetNPtBins(nPtBinLimits-1);
750 fnPtBinLimits = nPtBinLimits;
752 //cout<<"Changing also Global Index -> "<<fGlobalIndex<<endl;
753 fPtBinLimits = new Float_t[fnPtBinLimits];
754 for(Int_t ib=0; ib<nPtBinLimits; ib++) fPtBinLimits[ib]=ptBinLimits[ib];
758 //---------------------------------------------------------------------------
759 void AliRDHFCuts::SetVarNames(Int_t nVars,TString *varNames,Bool_t *isUpperCut){
760 // Set the variable names
765 //printf("Changing the variable names\n");
768 printf("Wrong number of variables: it has to be %d\n",fnVars);
772 fVarNames = new TString[nVars];
773 fIsUpperCut = new Bool_t[nVars];
774 for(Int_t iv=0; iv<nVars; iv++) {
775 fVarNames[iv] = varNames[iv];
776 fIsUpperCut[iv] = isUpperCut[iv];
781 //---------------------------------------------------------------------------
782 void AliRDHFCuts::SetVarsForOpt(Int_t nVars,Bool_t *forOpt) {
783 // Set the variables to be used for cuts optimization
786 delete [] fVarsForOpt;
788 //printf("Changing the variables for cut optimization\n");
791 if(nVars==0){//!=fnVars) {
792 printf("%d not accepted as number of variables: it has to be %d\n",nVars,fnVars);
797 fVarsForOpt = new Bool_t[fnVars];
798 for(Int_t iv=0; iv<fnVars; iv++) {
799 fVarsForOpt[iv]=forOpt[iv];
800 if(fVarsForOpt[iv]) fnVarsForOpt++;
806 //---------------------------------------------------------------------------
807 void AliRDHFCuts::SetUseCentrality(Int_t flag) {
809 // set centrality estimator
812 if(fUseCentrality<kCentOff||fUseCentrality>=kCentInvalid) AliWarning("Centrality estimator not valid");
818 //---------------------------------------------------------------------------
819 void AliRDHFCuts::SetCuts(Int_t nVars,Int_t nPtBins,Float_t **cutsRD) {
824 printf("Wrong number of variables: it has to be %d\n",fnVars);
827 if(nPtBins!=fnPtBins) {
828 printf("Wrong number of pt bins: it has to be %d\n",fnPtBins);
832 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
835 for(Int_t iv=0; iv<fnVars; iv++) {
837 for(Int_t ib=0; ib<fnPtBins; ib++) {
840 if(GetGlobalIndex(iv,ib)>=fGlobalIndex) {
841 cout<<"Overflow, exit..."<<endl;
845 fCutsRD[GetGlobalIndex(iv,ib)] = cutsRD[iv][ib];
851 //---------------------------------------------------------------------------
852 void AliRDHFCuts::SetCuts(Int_t glIndex,Float_t* cutsRDGlob){
856 if(glIndex != fGlobalIndex){
857 cout<<"Wrong array size: it has to be "<<fGlobalIndex<<endl;
860 if(!fCutsRD) fCutsRD = new Float_t[fGlobalIndex];
862 for(Int_t iGl=0;iGl<fGlobalIndex;iGl++){
863 fCutsRD[iGl] = cutsRDGlob[iGl];
867 //---------------------------------------------------------------------------
868 void AliRDHFCuts::PrintAll() const {
870 // print all cuts values
873 printf("Minimum vtx type %d\n",fMinVtxType);
874 printf("Minimum vtx contr %d\n",fMinVtxContr);
875 printf("Max vtx red chi2 %f\n",fMaxVtxRedChi2);
876 printf("Min SPD mult %d\n",fMinSPDMultiplicity);
877 printf("Use PID %d OldPid=%d\n",(Int_t)fUsePID,fPidHF ? fPidHF->GetOldPid() : -1);
878 printf("Remove daughters from vtx %d\n",(Int_t)fRemoveDaughtersFromPrimary);
879 printf("Recompute primary vertex %d\n",(Int_t)fRecomputePrimVertex);
880 printf("Physics selection: %s\n",fUsePhysicsSelection ? "Yes" : "No");
881 printf("Pileup rejection: %s\n",(fOptPileup > 0) ? "Yes" : "No");
882 if(fOptPileup==1) printf(" -- Reject pileup event");
883 if(fOptPileup==2) printf(" -- Reject tracks from pileup vtx");
884 if(fUseCentrality>0) {
885 TString estimator="";
886 if(fUseCentrality==1) estimator = "V0";
887 if(fUseCentrality==2) estimator = "Tracks";
888 if(fUseCentrality==3) estimator = "Tracklets";
889 if(fUseCentrality==4) estimator = "SPD clusters outer";
890 printf("Centrality class considered: %.1f-%.1f, estimated with %s",fMinCentrality,fMaxCentrality,estimator.Data());
892 if(fIsCandTrackSPDFirst) printf("Check for candidates with pt < %2.2f, that daughters fullfill kFirst criteria\n",fMaxPtCandTrackSPDFirst);
895 cout<<"Array of variables"<<endl;
896 for(Int_t iv=0;iv<fnVars;iv++){
897 cout<<fVarNames[iv]<<"\t";
902 cout<<"Array of optimization"<<endl;
903 for(Int_t iv=0;iv<fnVars;iv++){
904 cout<<fVarsForOpt[iv]<<"\t";
909 cout<<"Array of upper/lower cut"<<endl;
910 for(Int_t iv=0;iv<fnVars;iv++){
911 cout<<fIsUpperCut[iv]<<"\t";
916 cout<<"Array of ptbin limits"<<endl;
917 for(Int_t ib=0;ib<fnPtBinLimits;ib++){
918 cout<<fPtBinLimits[ib]<<"\t";
923 cout<<"Matrix of cuts"<<endl;
924 for(Int_t iv=0;iv<fnVars;iv++){
925 for(Int_t ib=0;ib<fnPtBins;ib++){
926 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<"\t";
935 //--------------------------------------------------------------------------
936 void AliRDHFCuts::PrintTrigger() const{
937 // print the trigger selection
939 printf("Selected trigger classes: %s %s\n",fTriggerClass[0].Data(),fTriggerClass[1].Data());
941 cout<<" Trigger selection pattern: ";
942 if( fTriggerMask & AliVEvent::kAny ) cout<<" kAny ";
943 if( fTriggerMask & AliVEvent::kAnyINT ) cout<<" kAnyINT ";
944 if( fTriggerMask & AliVEvent::kINT7 ) cout<<" kINT7 ";
945 if( fTriggerMask & AliVEvent::kMB ) cout<<" kMB ";
946 if( fTriggerMask & AliVEvent::kCINT5 ) cout<<" kCINT5 ";
947 if( fTriggerMask & AliVEvent::kCentral ) cout<<" kCentral ";
948 if( fTriggerMask & AliVEvent::kSemiCentral ) cout<<" kSemiCentral ";
949 if( fTriggerMask & AliVEvent::kEMCEGA ) cout<<" kEMCEGA ";
950 if( fTriggerMask & AliVEvent::kHighMult ) cout<<" kHighMult ";
951 if( fTriggerMask & AliVEvent::kFastOnly ) cout<<" kFastOnly ";
956 //---------------------------------------------------------------------------
957 void AliRDHFCuts::GetCuts(Float_t**& cutsRD) const{
962 //cout<<"Give back a "<<fnVars<<"x"<<fnPtBins<<" matrix."<<endl;
967 //cout<<"Initialization..."<<endl;
968 cutsRD=new Float_t*[fnVars];
969 for(iv=0; iv<fnVars; iv++) {
970 cutsRD[iv] = new Float_t[fnPtBins];
974 for(Int_t iGlobal=0; iGlobal<fGlobalIndex; iGlobal++) {
975 GetVarPtIndex(iGlobal,iv,ib);
976 cutsRD[iv][ib] = fCutsRD[iGlobal];
982 //---------------------------------------------------------------------------
983 Int_t AliRDHFCuts::GetGlobalIndex(Int_t iVar,Int_t iPtBin) const{
985 // give the global index from variable and pt bin
987 return iPtBin*fnVars+iVar;
990 //---------------------------------------------------------------------------
991 void AliRDHFCuts::GetVarPtIndex(Int_t iGlob, Int_t& iVar, Int_t& iPtBin) const {
993 //give the index of the variable and of the pt bin from the global index
995 iPtBin=(Int_t)iGlob/fnVars;
1001 //---------------------------------------------------------------------------
1002 Int_t AliRDHFCuts::PtBin(Double_t pt) const {
1004 //give the pt bin where the pt lies.
1007 if(pt<fPtBinLimits[0])return ptbin;
1008 for (Int_t i=0;i<fnPtBins;i++){
1009 if(pt<fPtBinLimits[i+1]) {
1016 //-------------------------------------------------------------------
1017 Float_t AliRDHFCuts::GetCutValue(Int_t iVar,Int_t iPtBin) const {
1019 // Give the value of cut set for the variable iVar and the pt bin iPtBin
1022 cout<<"Cuts not iniziaisez yet"<<endl;
1025 return fCutsRD[GetGlobalIndex(iVar,iPtBin)];
1027 //-------------------------------------------------------------------
1028 Float_t AliRDHFCuts::GetCentrality(AliAODEvent* aodEvent,AliRDHFCuts::ECentrality estimator) {
1030 // Get centrality percentile
1033 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aodEvent)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1034 if(mcArray) {fUseAOD049=kFALSE;}
1036 AliAODHeader *header=aodEvent->GetHeader();
1037 AliCentrality *centrality=header->GetCentralityP();
1039 Bool_t isSelRun=kFALSE;
1040 Int_t selRun[5]={138364, 138826, 138828, 138836, 138871};
1041 if(!centrality) return cent;
1043 if (estimator==kCentV0M){
1044 cent=(Float_t)(centrality->GetCentralityPercentile("V0M"));
1046 Int_t quality = centrality->GetQuality();
1047 if(quality<=1){ // fQuality==1 means rejected by zVertex cut that we apply a part and we want to keep separate (Giacomo)
1048 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1050 Int_t runnum=aodEvent->GetRunNumber();
1051 for(Int_t ir=0;ir<5;ir++){
1052 if(runnum==selRun[ir]){
1057 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
1061 //temporary fix for AOD049 outliers
1062 if(fUseAOD049&¢>=0){
1064 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1065 v0+=aodV0->GetMTotV0A();
1066 v0+=aodV0->GetMTotV0C();
1067 if(cent==0&&v0<19500)return -1;//filtering issue
1068 Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
1069 Float_t val= 1.30552 + 0.147931 * v0;
1070 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};
1071 if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )return -1;//outlier
1075 if (estimator==kCentTRK) {
1076 cent=(Float_t)(centrality->GetCentralityPercentile("TRK"));
1078 Int_t quality = centrality->GetQuality();
1080 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1082 Int_t runnum=aodEvent->GetRunNumber();
1083 for(Int_t ir=0;ir<5;ir++){
1084 if(runnum==selRun[ir]){
1089 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TRK");
1094 if (estimator==kCentTKL){
1095 cent=(Float_t)(centrality->GetCentralityPercentile("TKL"));
1097 Int_t quality = centrality->GetQuality();
1099 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1101 Int_t runnum=aodEvent->GetRunNumber();
1102 for(Int_t ir=0;ir<5;ir++){
1103 if(runnum==selRun[ir]){
1108 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("TKL");
1113 if (estimator==kCentCL1){
1114 cent=(Float_t)(centrality->GetCentralityPercentile("CL1"));
1116 Int_t quality = centrality->GetQuality();
1118 cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1120 Int_t runnum=aodEvent->GetRunNumber();
1121 for(Int_t ir=0;ir<5;ir++){
1122 if(runnum==selRun[ir]){
1127 if((quality==8||quality==9)&&isSelRun)cent=(Float_t)centrality->GetCentralityPercentileUnchecked("CL1");
1132 AliWarning("Centrality estimator not valid");
1141 //-------------------------------------------------------------------
1142 Bool_t AliRDHFCuts::CompareCuts(const AliRDHFCuts *obj) const {
1144 // Compare two cuts objects
1147 Bool_t areEqual=kTRUE;
1149 if(fMinVtxType!=obj->fMinVtxType) { printf("Minimum vtx type %d %d\n",fMinVtxType,obj->fMinVtxType); areEqual=kFALSE;}
1151 if(fMinVtxContr!=obj->fMinVtxContr) { printf("Minimum vtx contr %d %d\n",fMinVtxContr,obj->fMinVtxContr); areEqual=kFALSE;}
1153 if(TMath::Abs(fMaxVtxRedChi2-obj->fMaxVtxRedChi2)>1.e-10) { printf("Max vtx red chi2 %f %f\n",fMaxVtxRedChi2,obj->fMaxVtxRedChi2);areEqual=kFALSE;}
1155 if(fMinSPDMultiplicity!=obj->fMinSPDMultiplicity) { printf("Min SPD mult %d\n %d",fMinSPDMultiplicity,obj->fMinSPDMultiplicity);areEqual=kFALSE;}
1157 if(fUsePID!=obj->fUsePID) { printf("Use PID %d %d\n",(Int_t)fUsePID,(Int_t)obj->fUsePID); areEqual=kFALSE;}
1159 if(fRemoveDaughtersFromPrimary!=obj->fRemoveDaughtersFromPrimary) {printf("Remove daughters from vtx %d %d\n",(Int_t)fRemoveDaughtersFromPrimary,(Int_t)obj->fRemoveDaughtersFromPrimary); areEqual=kFALSE;}
1161 if(fTrackCuts->GetMinNClusterTPC()!=obj->fTrackCuts->GetMinNClusterTPC()) {printf("MinNClsTPC %d %d\n",fTrackCuts->GetMinNClusterTPC(),obj->fTrackCuts->GetMinNClusterTPC()); areEqual=kFALSE;}
1163 if(fTrackCuts->GetMinNClustersITS()!=obj->fTrackCuts->GetMinNClustersITS()) {printf("MinNClsITS %d %d\n",fTrackCuts->GetMinNClustersITS(),obj->fTrackCuts->GetMinNClustersITS()); areEqual=kFALSE;}
1165 if(TMath::Abs(fTrackCuts->GetMaxChi2PerClusterTPC()-obj->fTrackCuts->GetMaxChi2PerClusterTPC())>1.e-10) {printf("MaxChi2ClsTPC %f %f\n",fTrackCuts->GetMaxChi2PerClusterTPC(),obj->fTrackCuts->GetMaxChi2PerClusterTPC()); areEqual=kFALSE;}
1167 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;}
1171 for(Int_t iv=0;iv<fnVars;iv++) {
1172 for(Int_t ib=0;ib<fnPtBins;ib++) {
1173 if(TMath::Abs(fCutsRD[GetGlobalIndex(iv,ib)]-obj->fCutsRD[GetGlobalIndex(iv,ib)])>1.e-10) {
1174 cout<<"fCutsRD["<<iv<<"]["<<ib<<"] = "<<fCutsRD[GetGlobalIndex(iv,ib)]<<" "<<obj->fCutsRD[GetGlobalIndex(iv,ib)]<<"\n";
1183 //---------------------------------------------------------------------------
1184 void AliRDHFCuts::MakeTable() const {
1186 // print cuts values in table format
1189 TString ptString = "pT range";
1190 if(fVarNames && fPtBinLimits && fCutsRD){
1191 TString firstLine(Form("* %-15s",ptString.Data()));
1192 for (Int_t ivar=0; ivar<fnVars; ivar++){
1193 firstLine+=Form("* %-15s ",fVarNames[ivar].Data());
1194 if (ivar == fnVars){
1198 Printf("%s",firstLine.Data());
1200 for (Int_t ipt=0; ipt<fnPtBins; ipt++){
1202 if (ipt==fnPtBins-1){
1203 line=Form("* %5.1f < pt < inf ",fPtBinLimits[ipt]);
1206 line=Form("* %5.1f < pt < %4.1f ",fPtBinLimits[ipt],fPtBinLimits[ipt+1]);
1208 for (Int_t ivar=0; ivar<fnVars; ivar++){
1209 line+=Form("* %-15f ",fCutsRD[GetGlobalIndex(ivar,ipt)]);
1211 Printf("%s",line.Data());
1219 //--------------------------------------------------------------------------
1220 Bool_t AliRDHFCuts::RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d,
1221 AliAODEvent *aod) const
1224 // Recalculate primary vertex without daughters
1228 AliError("Can not remove daughters from vertex without AOD event");
1232 AliAODVertex *recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
1234 AliDebug(2,"Removal of daughter tracks failed");
1239 //set recalculed primary vertex
1240 d->SetOwnPrimaryVtx(recvtx);
1245 //--------------------------------------------------------------------------
1246 Bool_t AliRDHFCuts::SetMCPrimaryVtx(AliAODRecoDecayHF *d,AliAODEvent *aod) const
1249 // Recalculate primary vertex without daughters
1253 AliError("Can not get MC vertex without AOD event");
1258 AliAODMCHeader *mcHeader =
1259 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1261 AliError("Can not get MC vertex without AODMCHeader event");
1265 Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1266 mcHeader->GetVertex(pos);
1267 AliAODVertex *recvtx=new AliAODVertex(pos,covmatrix);
1270 AliDebug(2,"Removal of daughter tracks failed");
1274 //set recalculed primary vertex
1275 d->SetOwnPrimaryVtx(recvtx);
1277 d->RecalculateImpPars(recvtx,aod);
1283 //--------------------------------------------------------------------------
1284 void AliRDHFCuts::CleanOwnPrimaryVtx(AliAODRecoDecayHF *d,
1286 AliAODVertex *origownvtx) const
1289 // Clean-up own primary vertex if needed
1292 if(fRemoveDaughtersFromPrimary || fUseMCVertex) {
1293 d->UnsetOwnPrimaryVtx();
1295 d->SetOwnPrimaryVtx(origownvtx);
1296 delete origownvtx; origownvtx=NULL;
1298 d->RecalculateImpPars(d->GetPrimaryVtx(),aod);
1301 delete origownvtx; origownvtx=NULL;
1306 //--------------------------------------------------------------------------
1307 Bool_t AliRDHFCuts::IsSignalMC(AliAODRecoDecay *d,AliAODEvent *aod,Int_t pdg) const
1310 // Checks if this candidate is matched to MC signal
1313 if(!aod) return kFALSE;
1316 TClonesArray *mcArray = (TClonesArray*)((AliAODEvent*)aod)->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1318 if(!mcArray) return kFALSE;
1321 Int_t label = d->MatchToMC(pdg,mcArray);
1324 //printf("MATCH!\n");
1332 //--------------------------------------------------------------------------
1333 Bool_t AliRDHFCuts::RecomputePrimaryVertex(AliAODEvent* event) const{
1334 // recompute event primary vertex from AOD tracks
1336 AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());
1337 vertexer->SetITSMode();
1338 vertexer->SetMinClusters(3);
1340 AliAODVertex* pvtx=event->GetPrimaryVertex();
1341 if(strstr(pvtx->GetTitle(),"VertexerTracksWithConstraint")) {
1342 Float_t diamondcovxy[3];
1343 event->GetDiamondCovXY(diamondcovxy);
1344 Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
1345 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1346 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
1347 vertexer->SetVtxStart(diamond);
1348 delete diamond; diamond=NULL;
1351 AliESDVertex* vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
1352 if(!vertexESD) return kFALSE;
1353 if(vertexESD->GetNContributors()<=0) {
1354 //AliDebug(2,"vertexing failed");
1355 delete vertexESD; vertexESD=NULL;
1358 delete vertexer; vertexer=NULL;
1360 // convert to AliAODVertex
1361 Double_t pos[3],cov[6],chi2perNDF;
1362 vertexESD->GetXYZ(pos); // position
1363 vertexESD->GetCovMatrix(cov); //covariance matrix
1364 chi2perNDF = vertexESD->GetChi2toNDF();
1365 delete vertexESD; vertexESD=NULL;
1367 pvtx->SetPosition(pos[0],pos[1],pos[2]);
1368 pvtx->SetChi2perNDF(chi2perNDF);
1369 pvtx->SetCovMatrix(cov);