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