]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New methods for tagging decay channels
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliCFTaskVertexingHF.cxx
CommitLineData
379592af 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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//-----------------------------------------------------------------------
17// Class for HF corrections as a function of many variables
18// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19// Reco Acc + ITS Cl + PPR cuts
20// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
26//-----------------------------------------------------------------------
27//-----------------------------------------------------------------------
28// Base class for HF Unfolding (pt and eta)
29// correlation matrix filled at Acceptance and PPR level
30// Author: A.Grelli , Utrecht - agrelli@uu.nl
31//-----------------------------------------------------------------------
32#include <TCanvas.h>
33#include <TParticle.h>
34#include <TDatabasePDG.h>
f26726ed 35#include <TProfile.h>
379592af 36#include <TH1I.h>
37#include <TStyle.h>
38#include <TFile.h>
5b37c6e5 39#include <TF1.h>
379592af 40
41#include "AliCFTaskVertexingHF.h"
42#include "AliStack.h"
43#include "AliMCEvent.h"
44#include "AliCFManager.h"
45#include "AliCFContainer.h"
46#include "AliLog.h"
5cd139bc 47#include "AliInputEventHandler.h"
379592af 48#include "AliAnalysisManager.h"
49#include "AliAODHandler.h"
50#include "AliAODEvent.h"
51#include "AliAODRecoDecay.h"
52#include "AliAODRecoDecayHF.h"
53#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 54#include "AliAODRecoDecayHF3Prong.h"
55#include "AliAODRecoDecayHF4Prong.h"
56#include "AliAODRecoCascadeHF.h"
379592af 57#include "AliAODMCParticle.h"
58#include "AliAODMCHeader.h"
59#include "AliESDtrack.h"
60#include "TChain.h"
61#include "THnSparse.h"
62#include "TH2D.h"
63#include "AliESDtrackCuts.h"
f2703bd2 64#include "AliRDHFCuts.h"
379592af 65#include "AliRDHFCutsD0toKpi.h"
f2703bd2 66#include "AliRDHFCutsDplustoKpipi.h"
67#include "AliRDHFCutsDStartoKpipi.h"
68#include "AliRDHFCutsDstoKKpi.h"
69#include "AliRDHFCutsLctopKpi.h"
70#include "AliRDHFCutsD0toKpipipi.h"
5cd139bc 71#include "AliRDHFCutsLctoV0.h"
379592af 72#include "AliCFVertexingHF2Prong.h"
4e600f58 73#include "AliCFVertexingHF3Prong.h"
367e9aa3 74#include "AliCFVertexingHFCascade.h"
5cd139bc 75#include "AliCFVertexingHFLctoV0bachelor.h"
379592af 76#include "AliCFVertexingHF.h"
afb24c40 77#include "AliVertexingHFUtils.h"
f2703bd2 78#include "AliAnalysisDataSlot.h"
79#include "AliAnalysisDataContainer.h"
5cd139bc 80#include "AliPIDResponse.h"
379592af 81
82//__________________________________________________________________________
83AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
5cd139bc 84 AliAnalysisTaskSE(),
85 fCFManager(0x0),
86 fHistEventsProcessed(0x0),
87 fCorrelation(0x0),
f26726ed 88 fListProfiles(0),
5cd139bc 89 fCountMC(0),
90 fCountAcc(0),
91 fCountVertex(0),
92 fCountRefit(0),
93 fCountReco(0),
94 fCountRecoAcc(0),
95 fCountRecoITSClusters(0),
96 fCountRecoPPR(0),
97 fCountRecoPID(0),
98 fEvents(0),
99 fDecayChannel(0),
100 fFillFromGenerated(kFALSE),
101 fOriginDselection(0),
102 fAcceptanceUnf(kTRUE),
103 fCuts(0),
104 fUseWeight(kFALSE),
105 fWeight(1.),
106 fUseFlatPtWeight(kFALSE),
107 fUseZWeight(kFALSE),
108 fUseNchWeight(kFALSE),
06c14b1c 109 fUseTrackletsWeight(kFALSE),
5cd139bc 110 fNvar(0),
111 fPartName(""),
112 fDauNames(""),
113 fSign(2),
114 fCentralitySelection(kTRUE),
115 fFakeSelection(0),
116 fRejectIfNoQuark(kTRUE),
117 fUseMCVertex(kFALSE),
118 fDsOption(1),
119 fGenDsOption(3),
120 fConfiguration(kCheetah), // by default, setting the fast configuration
121 fFuncWeight(0x0),
122 fHistoMeasNch(0x0),
123 fHistoMCNch(0x0),
124 fResonantDecay(0),
125 fLctoV0bachelorOption(1),
ca7d039b 126 fGenLctoV0bachelorOption(0),
04b11915 127 fUseSelectionBit(kTRUE),
d74f7bc1 128 fPDGcode(0),
129 fMultiplicityEstimator(kNtrk10),
f26726ed 130 fRefMult(9.26),
131 fZvtxCorrectedNtrkEstimator(kFALSE),
d74f7bc1 132 fIsPPData(kFALSE)
379592af 133{
5cd139bc 134 //
135 //Default ctor
136 //
f26726ed 137 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
379592af 138}
139//___________________________________________________________________________
5b37c6e5 140AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
5cd139bc 141 AliAnalysisTaskSE(name),
142 fCFManager(0x0),
143 fHistEventsProcessed(0x0),
144 fCorrelation(0x0),
f26726ed 145 fListProfiles(0),
5cd139bc 146 fCountMC(0),
147 fCountAcc(0),
148 fCountVertex(0),
149 fCountRefit(0),
150 fCountReco(0),
151 fCountRecoAcc(0),
152 fCountRecoITSClusters(0),
153 fCountRecoPPR(0),
154 fCountRecoPID(0),
155 fEvents(0),
156 fDecayChannel(0),
157 fFillFromGenerated(kFALSE),
158 fOriginDselection(0),
159 fAcceptanceUnf(kTRUE),
160 fCuts(cuts),
161 fUseWeight(kFALSE),
162 fWeight(1.),
163 fUseFlatPtWeight(kFALSE),
164 fUseZWeight(kFALSE),
165 fUseNchWeight(kFALSE),
06c14b1c 166 fUseTrackletsWeight(kFALSE),
5cd139bc 167 fNvar(0),
168 fPartName(""),
169 fDauNames(""),
170 fSign(2),
171 fCentralitySelection(kTRUE),
172 fFakeSelection(0),
173 fRejectIfNoQuark(kTRUE),
174 fUseMCVertex(kFALSE),
175 fDsOption(1),
176 fGenDsOption(3),
177 fConfiguration(kCheetah), // by default, setting the fast configuration
178 fFuncWeight(func),
179 fHistoMeasNch(0x0),
180 fHistoMCNch(0x0),
181 fResonantDecay(0),
182 fLctoV0bachelorOption(1),
ca7d039b 183 fGenLctoV0bachelorOption(0),
04b11915 184 fUseSelectionBit(kTRUE),
d74f7bc1 185 fPDGcode(0),
186 fMultiplicityEstimator(kNtrk10),
f26726ed 187 fRefMult(9.26),
188 fZvtxCorrectedNtrkEstimator(kFALSE),
d74f7bc1 189 fIsPPData(kFALSE)
379592af 190{
5cd139bc 191 //
192 // Constructor. Initialization of Inputs and Outputs
193 //
194 /*
195 DefineInput(0) and DefineOutput(0)
196 are taken care of by AliAnalysisTaskSE constructor
197 */
198 DefineOutput(1,TH1I::Class());
199 DefineOutput(2,AliCFContainer::Class());
200 DefineOutput(3,THnSparseD::Class());
201 DefineOutput(4,AliRDHFCuts::Class());
f26726ed 202 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
203 DefineOutput(5,TList::Class()); // slot #5 keeps the zvtx Ntrakclets correction profiles
379592af 204
5cd139bc 205 fCuts->PrintAll();
379592af 206}
207
208//___________________________________________________________________________
209AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
210{
5cd139bc 211 //
212 // Assignment operator
213 //
214 if (this!=&c) {
215 AliAnalysisTaskSE::operator=(c) ;
216 fCFManager = c.fCFManager;
217 fHistEventsProcessed = c.fHistEventsProcessed;
218 fCuts = c.fCuts;
219 fFuncWeight = c.fFuncWeight;
220 fHistoMeasNch = c.fHistoMeasNch;
221 fHistoMCNch = c.fHistoMCNch;
f26726ed 222 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
5cd139bc 223 }
224 return *this;
379592af 225}
226
227//___________________________________________________________________________
228AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
5cd139bc 229 AliAnalysisTaskSE(c),
230 fCFManager(c.fCFManager),
231 fHistEventsProcessed(c.fHistEventsProcessed),
232 fCorrelation(c.fCorrelation),
f26726ed 233 fListProfiles(c.fListProfiles),
5cd139bc 234 fCountMC(c.fCountMC),
235 fCountAcc(c.fCountAcc),
236 fCountVertex(c.fCountVertex),
237 fCountRefit(c.fCountRefit),
238 fCountReco(c.fCountReco),
239 fCountRecoAcc(c.fCountRecoAcc),
240 fCountRecoITSClusters(c.fCountRecoITSClusters),
241 fCountRecoPPR(c.fCountRecoPPR),
242 fCountRecoPID(c.fCountRecoPID),
243 fEvents(c.fEvents),
244 fDecayChannel(c.fDecayChannel),
245 fFillFromGenerated(c.fFillFromGenerated),
246 fOriginDselection(c.fOriginDselection),
247 fAcceptanceUnf(c.fAcceptanceUnf),
248 fCuts(c.fCuts),
249 fUseWeight(c.fUseWeight),
250 fWeight(c.fWeight),
251 fUseFlatPtWeight(c.fUseFlatPtWeight),
252 fUseZWeight(c.fUseZWeight),
253 fUseNchWeight(c.fUseNchWeight),
06c14b1c 254 fUseTrackletsWeight(c.fUseTrackletsWeight),
5cd139bc 255 fNvar(c.fNvar),
256 fPartName(c.fPartName),
257 fDauNames(c.fDauNames),
258 fSign(c.fSign),
259 fCentralitySelection(c.fCentralitySelection),
260 fFakeSelection(c.fFakeSelection),
261 fRejectIfNoQuark(c.fRejectIfNoQuark),
262 fUseMCVertex(c.fUseMCVertex),
263 fDsOption(c.fDsOption),
264 fGenDsOption(c.fGenDsOption),
265 fConfiguration(c.fConfiguration),
266 fFuncWeight(c.fFuncWeight),
267 fHistoMeasNch(c.fHistoMeasNch),
268 fHistoMCNch(c.fHistoMCNch),
269 fResonantDecay(c.fResonantDecay),
270 fLctoV0bachelorOption(c.fLctoV0bachelorOption),
ca7d039b 271 fGenLctoV0bachelorOption(c.fGenLctoV0bachelorOption),
04b11915 272 fUseSelectionBit(c.fUseSelectionBit),
d74f7bc1 273 fPDGcode(c.fPDGcode),
274 fMultiplicityEstimator(c.fMultiplicityEstimator),
f26726ed 275 fRefMult(c.fRefMult),
276 fZvtxCorrectedNtrkEstimator(c.fZvtxCorrectedNtrkEstimator),
d74f7bc1 277 fIsPPData(c.fIsPPData)
379592af 278{
5cd139bc 279 //
280 // Copy Constructor
281 //
f26726ed 282 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
379592af 283}
284
285//___________________________________________________________________________
f2703bd2 286AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
287{
5cd139bc 288 //
289 //destructor
290 //
291 if (fCFManager) delete fCFManager ;
292 if (fHistEventsProcessed) delete fHistEventsProcessed ;
293 if (fCorrelation) delete fCorrelation ;
f26726ed 294 if (fListProfiles) delete fListProfiles;
5cd139bc 295 if (fCuts) delete fCuts;
296 if (fFuncWeight) delete fFuncWeight;
297 if (fHistoMeasNch) delete fHistoMeasNch;
298 if (fHistoMCNch) delete fHistoMCNch;
f26726ed 299 for(Int_t i=0; i<4; i++) { if(fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i]; }
379592af 300}
301
f2703bd2 302//_________________________________________________________________________-
303void AliCFTaskVertexingHF::Init()
304{
5cd139bc 305 //
306 // Initialization
307 //
f2703bd2 308
5cd139bc 309 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
310 if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
f26726ed 311 if(fUseWeight && fUseNchWeight) { AliInfo("Beware, using at the same time pt and Nch weights, please check"); }
5cd139bc 312 if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; }
313 if(fUseNchWeight) CreateMeasuredNchHisto();
8a983038 314
5cd139bc 315 AliRDHFCuts *copyfCuts = 0x0;
316 if (!fCuts){
317 AliFatal("No cuts defined - Exiting...");
318 return;
319 }
320
321 switch (fDecayChannel){
322 case 2:{
04b11915 323 fPDGcode = 421;
5cd139bc 324 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
325 switch (fConfiguration) {
326 case kSnail: // slow configuration: all variables in
327 fNvar = 16;
328 break;
329 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
330 fNvar = 8;
331 break;
332 }
333 fPartName="D0";
334 fDauNames="K+pi";
335 break;
336 }
337 case 21:{
04b11915 338 fPDGcode = 413;
5cd139bc 339 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
340 switch (fConfiguration) {
341 case kSnail: // slow configuration: all variables in
342 fNvar = 16;
343 break;
344 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
345 fNvar = 8;
346 break;
347 }
348 fPartName="Dstar";
349 fDauNames="K+pi+pi";
350 break;
351 }
352 case 22:{
04b11915 353 fPDGcode = 4122;
5cd139bc 354 copyfCuts = new AliRDHFCutsLctoV0(*(static_cast<AliRDHFCutsLctoV0*>(fCuts)));
355 switch (fConfiguration) {
356 case kSnail: // slow configuration: all variables in
357 fNvar = 16;
358 break;
359 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
360 fNvar = 8;
361 break;
362 }
363 fPartName="Lambdac";
364 fDauNames="V0+bachelor";
365 break;
366 }
367 case 31:{
04b11915 368 fPDGcode = 411;
5cd139bc 369 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
370 switch (fConfiguration) {
371 case kSnail: // slow configuration: all variables in
372 fNvar = 14;
373 break;
374 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
375 fNvar = 8;
376 break;
377 }
378 fPartName="Dplus";
379 fDauNames="K+pi+pi";
380 break;
381 }
382 case 32:{
04b11915 383 fPDGcode = 4122;
5cd139bc 384 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
385 switch (fConfiguration) {
386 case kSnail: // slow configuration: all variables in
387 fNvar = 18;
388 break;
389 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
390 fNvar = 8;
391 break;
392 }
393 fPartName="Lambdac";
394 fDauNames="p+K+pi";
395 break;
396 }
397 case 33:{
04b11915 398 fPDGcode = 431;
5cd139bc 399 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
400 switch (fConfiguration) {
401 case kSnail: // slow configuration: all variables in
402 fNvar = 14;
403 break;
404 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
405 fNvar = 8;
406 break;
407 }
408 fPartName="Ds";
409 fDauNames="K+K+pi";
410 break;
411 }
412 case 4:{
04b11915 413 fPDGcode = 421;
5cd139bc 414 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
415 switch (fConfiguration) {
416 case kSnail: // slow configuration: all variables in
417 fNvar = 16;
418 break;
419 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
420 fNvar = 8;
421 break;
422 }
423 fPartName="D0";
424 fDauNames="K+pi+pi+pi";
425 break;
426 }
427 default:
428 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
429 break;
430 }
f2703bd2 431
5cd139bc 432 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
433 if (copyfCuts){
434 copyfCuts->SetName(nameoutput);
8a983038 435
5cd139bc 436 //Post the data
437 PostData(4, copyfCuts);
438 }
439 else{
440 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
f26726ed 441 }
442
443 fListProfiles = new TList();
444 fListProfiles->SetOwner();
445 TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
446 for(Int_t i=0; i<4; i++){
447 if(fMultEstimatorAvg[i]){
448 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
449 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
450 fListProfiles->Add(hprof);
451 }
452 }
453 PostData(5,fListProfiles);
f2703bd2 454
5cd139bc 455 return;
f2703bd2 456}
457
379592af 458//_________________________________________________
459void AliCFTaskVertexingHF::UserExec(Option_t *)
460{
5cd139bc 461 //
462 // Main loop function
463 //
379592af 464
5cd139bc 465 PostData(1,fHistEventsProcessed) ;
466 PostData(2,fCFManager->GetParticleContainer()) ;
467 PostData(3,fCorrelation) ;
b7af2639 468
5cd139bc 469 if (fFillFromGenerated){
470 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
471 }
379592af 472
5cd139bc 473 if (!fInputEvent) {
474 Error("UserExec","NO EVENT FOUND!");
475 return;
476 }
379592af 477
5cd139bc 478 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
379592af 479
5cd139bc 480 TClonesArray *arrayBranch=0;
379592af 481
5cd139bc 482 if(!aodEvent && AODEvent() && IsStandardAOD()) {
483 // In case there is an AOD handler writing a standard AOD, use the AOD
484 // event in memory rather than the input (ESD) event.
485 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
486 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
487 // have to taken from the AOD event hold by the AliAODExtension
488 AliAODHandler* aodHandler = (AliAODHandler*)
489 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
490 if(aodHandler->GetExtensions()) {
491 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
492 AliAODEvent *aodFromExt = ext->GetAOD();
f2703bd2 493
5cd139bc 494 switch (fDecayChannel){
495 case 2:{
496 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
497 break;
498 }
499 case 21:{
500 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
501 break;
502 }
503 case 22:{
504 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
505 break;
506 }
507 case 31:
508 case 32:
509 case 33:{
510 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
511 break;
512 }
513 case 4:{
514 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
515 break;
516 }
517 default:
518 break;
519 }
520 }
521 }
522 else {
523 switch (fDecayChannel){
524 case 2:{
525 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
526 break;
527 }
528 case 21:{
529 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
530 break;
531 }
532 case 22:{
533 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
534 break;
535 }
536 case 31:
537 case 32:
538 case 33:{
539 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
540 break;
541 }
542 case 4:{
543 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
544 break;
545 }
546 default:
547 break;
548 }
549 }
ec5288c3 550
5cd139bc 551 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
552 if (!aodVtx) return;
f2703bd2 553
5cd139bc 554 if (!arrayBranch) {
555 AliError("Could not find array of HF vertices");
556 return;
557 }
379592af 558
5cd139bc 559 fEvents++;
6e2e4f50 560
5cd139bc 561 fCFManager->SetRecEventInfo(aodEvent);
562 fCFManager->SetMCEventInfo(aodEvent);
379592af 563
5cd139bc 564 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 565
5cd139bc 566 //loop on the MC event
379592af 567
5cd139bc 568 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
569 if (!mcArray) {
570 AliError("Could not find Monte-Carlo in AOD");
571 return;
572 }
573 Int_t icountMC = 0;
574 Int_t icountAcc = 0;
575 Int_t icountReco = 0;
576 Int_t icountVertex = 0;
577 Int_t icountRefit = 0;
578 Int_t icountRecoAcc = 0;
579 Int_t icountRecoITSClusters = 0;
580 Int_t icountRecoPPR = 0;
581 Int_t icountRecoPID = 0;
582 Int_t cquarks = 0;
f2703bd2 583
5cd139bc 584 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
585 if (!mcHeader) {
586 AliError("Could not find MC Header in AOD");
587 return;
588 }
e11ae259 589
78471b33 590 fHistEventsProcessed->Fill(0.5);
591
5cd139bc 592 Double_t* containerInput = new Double_t[fNvar];
593 Double_t* containerInputMC = new Double_t[fNvar];
e11ae259 594
f2703bd2 595
5cd139bc 596 AliCFVertexingHF* cfVtxHF=0x0;
597 switch (fDecayChannel){
598 case 2:{
599 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
600 break;
601 }
602 case 21:{
603 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
604 break;
605 }
606 case 22:{
607 cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption); // Lc -> K0S+proton
608 break;
609 }
610 case 31:
611 // case 32:
612 case 33:{
613 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
614 if(fDecayChannel==33){
615 ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption);
616 }
617 break;
618 }
619 case 32:{
620 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel,fResonantDecay);
621 }
622 case 4:{
623 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
624 break;
625 }
626 default:
627 break;
628 }
629 if (!cfVtxHF){
630 AliError("No AliCFVertexingHF initialized");
631 delete[] containerInput;
632 delete[] containerInputMC;
633 return;
634 }
f2703bd2 635
5cd139bc 636 Double_t zPrimVertex = aodVtx ->GetZ();
637 Double_t zMCVertex = mcHeader->GetVtxZ();
638 Int_t runnumber = aodEvent->GetRunNumber();
f26726ed 639
640 // Multiplicity definition with tracklets
641 Double_t nTracklets = 0;
642 Int_t nTrackletsEta10 = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.);
643 Int_t nTrackletsEta16 = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.6,1.6);
644 nTracklets = (Double_t)nTrackletsEta10;
645 if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16 - nTrackletsEta10); }
646
647 // Apply the Ntracklets z-vtx data driven correction
648 if(fZvtxCorrectedNtrkEstimator) {
649 TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
650 if(estimatorAvg) {
651 Int_t nTrackletsEta10Corr = AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta10,zPrimVertex,fRefMult);
652 Int_t nTrackletsEta16Corr = AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta16,zPrimVertex,fRefMult);
653 nTracklets = (Double_t)nTrackletsEta10Corr;
654 if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16Corr - nTrackletsEta10Corr); }
655 }
656 }
657
658
5cd139bc 659 fWeight=1.;
660 if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber);
661 if(fUseNchWeight){
662 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
4fa6d5f8 663 if(!fUseTrackletsWeight) fWeight *= GetNchWeight(nChargedMCPhysicalPrimary);
664 else fWeight *= GetNchWeight(nTracklets);
665 AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight));
5cd139bc 666 }
f26726ed 667 Double_t eventWeight=fWeight;
17bf1a34 668
5cd139bc 669 if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
670 AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
671 delete[] containerInput;
672 delete[] containerInputMC;
673 return;
674 }
751b8274 675
8bbb4282 676 if(aodEvent->GetTriggerMask()==0 &&
677 (runnumber>=195344 && runnumber<=195677)){
678 AliDebug(3,"Event rejected because of null trigger mask");
679 delete[] containerInput;
680 delete[] containerInputMC;
681 return;
682 }
683
5cd139bc 684 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
685 if (fDecayChannel == 21){
686 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
687 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
688 trackCuts[iProng]=fCuts->GetTrackCuts();
689 }
690 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
691 }
692 else if (fDecayChannel == 22) {
693 // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters
694 trackCuts[0]=fCuts->GetTrackCuts();
695 trackCuts[1]=fCuts->GetTrackCutsV0daughters();
696 trackCuts[2]=fCuts->GetTrackCutsV0daughters();
697 }
698 else {
699 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
700 trackCuts[iProng]=fCuts->GetTrackCuts();
701 }
702 }
2bf2e62b 703
5cd139bc 704 //General settings: vertex, feed down and fill reco container with generated values.
705 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
706 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
707 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
708 cfVtxHF->SetNVar(fNvar);
709 cfVtxHF->SetFakeSelection(fFakeSelection);
710 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
711 cfVtxHF->SetConfiguration(fConfiguration);
712
713 // switch-off the trigger class selection (doesn't work for MC)
714 fCuts->SetTriggerClass("");
715
716 // MC vertex, to be used, in case, for pp
717 if (fUseMCVertex) fCuts->SetUseMCVertex();
718
719 if (fCentralitySelection){ // keep only the requested centrality
720 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
721 delete[] containerInput;
722 delete[] containerInputMC;
723 delete [] trackCuts;
724 return;
725 }
726 } else { // keep all centralities
727 fCuts->SetMinCentrality(0.);
728 fCuts->SetMaxCentrality(100.);
729 }
8b1b055d 730
d74f7bc1 731 Float_t centValue = 0.;
732 if(!fIsPPData) centValue = fCuts->GetCentrality(aodEvent);
5cd139bc 733 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 734
f26726ed 735 // multiplicity estimator with VZERO
d74f7bc1 736 Double_t vzeroMult=0;
737 AliAODVZERO *vzeroAOD = (AliAODVZERO*)aodEvent->GetVZEROData();
738 if(vzeroAOD) vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
739
f26726ed 740 Double_t multiplicity = nTracklets; // set to the Ntracklet estimator
d74f7bc1 741 if(fMultiplicityEstimator==kVZERO) { multiplicity = vzeroMult; }
742
743
5cd139bc 744 cfVtxHF->SetMultiplicity(multiplicity);
d74f7bc1 745
746 // printf("Multiplicity estimator %d, value %2.2f\n",fMultiplicityEstimator,multiplicity);
ec5288c3 747
5cd139bc 748 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
749 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
750 if (!mcPart){
751 AliError("Failed casting particle from MC array!, Skipping particle");
752 continue;
753 }
04b11915 754
755 //counting c quarks
756 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
757
5cd139bc 758 // check the MC-level cuts, must be the desidered particle
759 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
760 AliDebug(2,"Check the MC-level cuts - not desidered particle");
761 continue; // 0 stands for MC level
762 }
763 cfVtxHF->SetMCCandidateParam(iPart);
b7af2639 764
b7af2639 765
5cd139bc 766 if (!(cfVtxHF->SetLabelArray())){
767 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
768 continue;
769 }
2bf2e62b 770
5cd139bc 771 //check the candiate family at MC level
772 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
773 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
774 continue;
775 }
776 else{
d74f7bc1 777 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
5cd139bc 778 }
f2703bd2 779
5cd139bc 780 //Fill the MC container
781 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
782 AliDebug(2,Form("mcContainerFilled = %d)",mcContainerFilled));
783 if (mcContainerFilled) {
784 if (fUseWeight){
785 if (fFuncWeight){ // using user-defined function
786 AliDebug(2,"Using function");
f26726ed 787 fWeight = eventWeight*fFuncWeight->Eval(containerInputMC[0]);
5cd139bc 788 }
789 else{ // using FONLL
790 AliDebug(2,"Using FONLL");
f26726ed 791 fWeight = eventWeight*GetWeight(containerInputMC[0]);
5cd139bc 792 }
793 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
794 }
795 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
796 //MC Limited Acceptance
797 if (TMath::Abs(containerInputMC[1]) < 0.5) {
798 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
799 AliDebug(3,"MC Lim Acc container filled\n");
800 }
f2703bd2 801
5cd139bc 802 //MC
803 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
804 icountMC++;
805 AliDebug(3,"MC cointainer filled \n");
f2703bd2 806
5cd139bc 807 // MC in acceptance
808 // check the MC-Acceptance level cuts
809 // since standard CF functions are not applicable, using Kine Cuts on daughters
810 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
811 if (mcAccepStep){
812 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
813 AliDebug(3,"MC acceptance cut passed\n");
814 icountAcc++;
f2703bd2 815
5cd139bc 816 //MC Vertex step
817 if (fCuts->IsEventSelected(aodEvent)){
818 // filling the container if the vertex is ok
819 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
820 AliDebug(3,"Vertex cut passed and container filled\n");
821 icountVertex++;
f2703bd2 822
5cd139bc 823 //mc Refit requirement
824 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
825 if (mcRefitStep){
826 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
827 AliDebug(3,"MC Refit cut passed and container filled\n");
828 icountRefit++;
829 }
830 else{
831 AliDebug(3,"MC Refit cut not passed\n");
832 continue;
833 }
379592af 834 }
5cd139bc 835 else{
836 AliDebug (3, "MC vertex step not passed\n");
837 continue;
838 }
839 }
840 else{
841 AliDebug (3, "MC in acceptance step not passed\n");
842 continue;
843 }
844 }
845 else {
846 AliDebug (3, "MC container not filled\n");
847 }
848 }
f2703bd2 849
5cd139bc 850 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
851 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
852 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
853 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
854 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
855
856 // Now go to rec level
857 fCountMC += icountMC;
858 fCountAcc += icountAcc;
859 fCountVertex+= icountVertex;
860 fCountRefit+= icountRefit;
861
862 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 863
5cd139bc 864 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
865 AliAODRecoDecayHF* charmCandidate=0x0;
866 switch (fDecayChannel){
867 case 2:{
868 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
869 break;
870 }
871 case 21:
872 case 22:{
873 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
874 break;
875 }
876 case 31:
877 case 32:
878 case 33:{
879 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
880 break;
881 }
882 case 4:{
883 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
884 break;
885 }
886 default:
887 break;
888 }
f2703bd2 889
5cd139bc 890 Bool_t unsetvtx=kFALSE;
891 if(!charmCandidate->GetOwnPrimaryVtx()) {
892 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
893 unsetvtx=kTRUE;
894 }
f2703bd2 895
5cd139bc 896 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
897 if (!signAssociation){
898 charmCandidate = 0x0;
899 continue;
900 }
3ee5eb83 901
5cd139bc 902 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
903 if (isPartOrAntipart == 0){
904 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
905 continue;
906 }
3ee5eb83 907
5cd139bc 908 AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
6694a2a8 909
5cd139bc 910 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
911 if (recoContFilled){
6694a2a8 912
5cd139bc 913 // weight according to pt
914 if (fUseWeight){
915 if (fFuncWeight){ // using user-defined function
916 AliDebug(2, "Using function");
f26726ed 917 fWeight = eventWeight*fFuncWeight->Eval(containerInput[0]);
5cd139bc 918 }
919 else{ // using FONLL
920 AliDebug(2, "Using FONLL");
f26726ed 921 fWeight = eventWeight*GetWeight(containerInput[0]);
5cd139bc 922 }
923 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
924 }
6694a2a8 925
5cd139bc 926 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 927
5cd139bc 928 //Reco Step
929 Bool_t recoStep = cfVtxHF->RecoStep();
930 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 931
5cd139bc 932
ca7d039b 933 // Selection on the filtering bit
934 Bool_t isBitSelected = kTRUE;
935 if(fDecayChannel==2) {
936 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
937 }else if(fDecayChannel==31){
938 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
939 }else if(fDecayChannel==33){
940 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
941 }
942 if(!isBitSelected) continue;
943
944
945
5cd139bc 946 if (recoStep && recoContFilled && vtxCheck){
947 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
948 icountReco++;
949 AliDebug(3,"Reco step passed and container filled\n");
f2703bd2 950
5cd139bc 951 //Reco in the acceptance -- take care of UNFOLDING!!!!
952 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
953 if (recoAcceptanceStep) {
954 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
955 icountRecoAcc++;
956 AliDebug(3,"Reco acceptance cut passed and container filled\n");
f2703bd2 957
5cd139bc 958 if(fAcceptanceUnf){
959 Double_t fill[4]; //fill response matrix
04b11915 960 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 961 if (bUnfolding) fCorrelation->Fill(fill);
962 }
f2703bd2 963
5cd139bc 964 //Number of ITS cluster requirements
965 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
966 if (recoITSnCluster){
967 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
968 icountRecoITSClusters++;
969 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
f2703bd2 970
5cd139bc 971 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
972 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
973 fCuts->SetUsePID(kFALSE);
974 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
975
976 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
977 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
978 if(keepDs) recoAnalysisCuts=3;
979 }
980 else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
981 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
982 if (keepLctoV0bachelor) recoAnalysisCuts=3;
983 }
984
6c62cb59 985
379592af 986
5cd139bc 987 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
988 Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
989 if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
fdd5a897 990
5cd139bc 991 if (tempAn){
992 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
993 icountRecoPPR++;
994 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
995 //pid selection
996 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
997 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
998 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
999
1000 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1001 Bool_t keepDs=ProcessDs(recoPidSelection);
1002 if(keepDs) recoPidSelection=3;
1003 } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
04b11915 1004 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1005 if (keepLctoV0bachelor) recoPidSelection=3;
5cd139bc 1006 }
1007
1008 Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1009 if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1010
1011 if (tempPid){
1012 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
1013 icountRecoPID++;
1014 AliDebug(3,"Reco PID cuts passed and container filled \n");
1015 if(!fAcceptanceUnf){
1016 Double_t fill[4]; //fill response matrix
04b11915 1017 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 1018 if (bUnfolding) fCorrelation->Fill(fill);
379592af 1019 }
5cd139bc 1020 }
1021 else {
1022 AliDebug(3, "Analysis Cuts step not passed \n");
1023 continue;
1024 }
1025 }
1026 else {
1027 AliDebug(3, "PID selection not passed \n");
1028 continue;
1029 }
1030 }
1031 else{
1032 AliDebug(3, "Number of ITS cluster step not passed\n");
1033 continue;
1034 }
1035 }
1036 else{
1037 AliDebug(3, "Reco acceptance step not passed\n");
1038 continue;
1039 }
1040 }
1041 else {
1042 AliDebug(3, "Reco step not passed\n");
1043 continue;
1044 }
1045 }
379592af 1046
5cd139bc 1047 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1048 } // end loop on candidate
f2703bd2 1049
5cd139bc 1050 fCountReco+= icountReco;
1051 fCountRecoAcc+= icountRecoAcc;
1052 fCountRecoITSClusters+= icountRecoITSClusters;
1053 fCountRecoPPR+= icountRecoPPR;
1054 fCountRecoPID+= icountRecoPID;
379592af 1055
78471b33 1056 fHistEventsProcessed->Fill(1.5);
5cd139bc 1057
1058 delete[] containerInput;
1059 delete[] containerInputMC;
1060 delete cfVtxHF;
1061 if (trackCuts){
1062 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1063 // delete [] trackCuts[i];
1064 // }
1065 delete [] trackCuts;
1066 }
ec5288c3 1067
1068
379592af 1069}
1070
1071//___________________________________________________________________________
1072void AliCFTaskVertexingHF::Terminate(Option_t*)
1073{
5cd139bc 1074 // The Terminate() function is the last function to be called during
1075 // a query. It always runs on the client, it can be used to present
1076 // the results graphically or save the results to file.
379592af 1077
5cd139bc 1078 AliAnalysisTaskSE::Terminate();
1079
1080 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1081 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1082 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
1083 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
1084 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1085 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
1086 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
1087 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
1088 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
379592af 1089
5cd139bc 1090 // draw some example plots....
1091 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1092 if(!cont) {
1093 printf("CONTAINER NOT FOUND\n");
1094 return;
1095 }
1096 // projecting the containers to obtain histograms
1097 // first argument = variable, second argument = step
ec5288c3 1098
5cd139bc 1099 TH1D** h = new TH1D*[3];
1100 Int_t nvarToPlot = 0;
1101 if (fConfiguration == kSnail){
1102 //h = new TH1D[3][12];
1103 for (Int_t ih = 0; ih<3; ih++){
1104 if(fDecayChannel==22){
1105 nvarToPlot = 16;
1106 } else {
1107 nvarToPlot = 12;
1108 }
1109 h[ih] = new TH1D[nvarToPlot];
1110 }
1111 for(Int_t iC=1;iC<nvarToPlot; iC++){
1112 // MC-level
1113 h[0][iC] = *(cont->ShowProjection(iC,0));
1114 // MC-Acceptance level
1115 h[1][iC] = *(cont->ShowProjection(iC,1));
1116 // Reco-level
1117 h[2][iC] = *(cont->ShowProjection(iC,4));
1118 }
1119 }
1120 else {
1121 //h = new TH1D[3][12];
1122 nvarToPlot = 8;
1123 for (Int_t ih = 0; ih<3; ih++){
1124 h[ih] = new TH1D[nvarToPlot];
1125 }
1126 for(Int_t iC=0;iC<nvarToPlot; iC++){
1127 // MC-level
1128 h[0][iC] = *(cont->ShowProjection(iC,0));
1129 // MC-Acceptance level
1130 h[1][iC] = *(cont->ShowProjection(iC,1));
1131 // Reco-level
1132 h[2][iC] = *(cont->ShowProjection(iC,4));
1133 }
1134 }
1135 TString* titles;
1136 //Int_t nvarToPlot = 0;
1137 if (fConfiguration == kSnail){
1138 if(fDecayChannel==31){
1139 //nvarToPlot = 12;
1140 titles = new TString[nvarToPlot];
1141 titles[0]="pT_Dplus (GeV/c)";
1142 titles[1]="rapidity";
1143 titles[2]="phi (rad)";
1144 titles[3]="cT (#mum)";
1145 titles[4]="cosPointingAngle";
1146 titles[5]="pT_1 (GeV/c)";
1147 titles[6]="pT_2 (GeV/c)";
1148 titles[7]="pT_3 (GeV/c)";
1149 titles[8]="d0_1 (#mum)";
1150 titles[9]="d0_2 (#mum)";
1151 titles[10]="d0_3 (#mum)";
1152 titles[11]="zVertex (cm)";
1153 } else if (fDecayChannel==22) {
1154 //nvarToPlot = 16;
1155 titles = new TString[nvarToPlot];
1156 titles[0]="pT_Lc (GeV/c)";
1157 titles[1]="rapidity";
1158 titles[2]="phi (rad)";
1159 titles[3]="cosPAV0";
1160 titles[4]="onTheFlyStatusV0";
1161 titles[5]="centrality";
1162 titles[6]="fake";
1163 titles[7]="multiplicity";
1164 titles[8]="pT_bachelor (GeV/c)";
1165 titles[9]="pT_V0pos (GeV/c)";
1166 titles[10]="pT_V0neg (GeV/c)";
1167 titles[11]="invMassV0 (GeV/c2)";
1168 titles[12]="dcaV0 (nSigma)";
1169 titles[13]="c#tauV0 (#mum)";
1170 titles[14]="c#tau (#mum)";
1171 titles[15]="cosPA";
1172 } else {
1173 //nvarToPlot = 12;
1174 titles = new TString[nvarToPlot];
1175 titles[0]="pT_D0 (GeV/c)";
1176 titles[1]="rapidity";
1177 titles[2]="cosThetaStar";
1178 titles[3]="pT_pi (GeV/c)";
1179 titles[4]="pT_K (Gev/c)";
1180 titles[5]="cT (#mum)";
1181 titles[6]="dca (#mum)";
1182 titles[7]="d0_pi (#mum)";
1183 titles[8]="d0_K (#mum)";
1184 titles[9]="d0xd0 (#mum^2)";
1185 titles[10]="cosPointingAngle";
1186 titles[11]="phi (rad)";
1187 }
1188 }
1189 else {
1190 //nvarToPlot = 8;
1191 titles = new TString[nvarToPlot];
1192 if (fDecayChannel==22) {
1193 titles[0]="pT_candidate (GeV/c)";
1194 titles[1]="rapidity";
1195 titles[2]="phi (rad)";
1196 titles[3]="cosPAV0";
1197 titles[4]="onTheFlyStatusV0";
1198 titles[5]="centrality";
1199 titles[6]="fake";
1200 titles[7]="multiplicity";
1201 } else {
1202 titles[0]="pT_candidate (GeV/c)";
1203 titles[1]="rapidity";
1204 titles[2]="cT (#mum)";
1205 titles[3]="phi";
1206 titles[4]="z_{vtx}";
1207 titles[5]="centrality";
1208 titles[6]="fake";
1209 titles[7]="multiplicity";
1210 }
1211 }
1212
1213 Int_t markers[16]={20,24,21,25,27,28,
1214 20,24,21,25,27,28,
1215 20,24,21,25};
1216 Int_t colors[3]={2,8,4};
1217 for(Int_t iC=0;iC<nvarToPlot; iC++){
1218 for(Int_t iStep=0;iStep<3;iStep++){
1219 h[iStep][iC].SetTitle(titles[iC].Data());
1220 h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1221 Double_t maxh=h[iStep][iC].GetMaximum();
1222 h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1223 h[iStep][iC].SetMarkerStyle(markers[iC]);
1224 h[iStep][iC].SetMarkerColor(colors[iStep]);
1225 }
1226 }
379592af 1227
5cd139bc 1228 gStyle->SetCanvasColor(0);
1229 gStyle->SetFrameFillColor(0);
1230 gStyle->SetTitleFillColor(0);
1231 gStyle->SetStatColor(0);
379592af 1232
5cd139bc 1233 // drawing in 2 separate canvas for a matter of clearity
1234 TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1235 c1->Divide(3,4);
1236 Int_t iPad=1;
1237 for(Int_t iVar=0; iVar<4; iVar++){
1238 c1->cd(iPad++);
1239 h[0][iVar].DrawCopy("p");
1240 c1->cd(iPad++);
1241 h[1][iVar].DrawCopy("p");
1242 c1->cd(iPad++);
1243 h[2][iVar].DrawCopy("p");
1244 }
379592af 1245
5cd139bc 1246 TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1247 c2->Divide(3,4);
1248 iPad=1;
1249 for(Int_t iVar=4; iVar<8; iVar++){
1250 c2->cd(iPad++);
1251 h[0][iVar].DrawCopy("p");
1252 c2->cd(iPad++);
1253 h[1][iVar].DrawCopy("p");
1254 c2->cd(iPad++);
1255 h[2][iVar].DrawCopy("p");
1256 }
31da6b05 1257
5cd139bc 1258 if (fConfiguration == kSnail){
1259 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1260 c3->Divide(3,4);
1261 iPad=1;
1262 for(Int_t iVar=8; iVar<12; iVar++){
1263 c3->cd(iPad++);
1264 h[0][iVar].DrawCopy("p");
1265 c3->cd(iPad++);
1266 h[1][iVar].DrawCopy("p");
1267 c3->cd(iPad++);
1268 h[2][iVar].DrawCopy("p");
1269 }
1270 if (fDecayChannel==22) {
1271 TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1272 c4->Divide(3,4);
1273 iPad=1;
1274 for(Int_t iVar=12; iVar<16; iVar++){
1275 c4->cd(iPad++);
1276 h[0][iVar].DrawCopy("p");
1277 c4->cd(iPad++);
1278 h[1][iVar].DrawCopy("p");
1279 c4->cd(iPad++);
1280 h[2][iVar].DrawCopy("p");
1281 }
1282 }
1283 }
31da6b05 1284
5cd139bc 1285 /*
1286 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
379592af 1287
5cd139bc 1288 TH2D* corr1 = hcorr->Projection(0,2);
1289 TH2D* corr2 = hcorr->Projection(1,3);
379592af 1290
5cd139bc 1291 TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1292 c7->Divide(2,1);
1293 c7->cd(1);
1294 corr1->DrawCopy("text");
1295 c7->cd(2);
1296 corr2->DrawCopy("text");
1297 */
1298 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
379592af 1299
5cd139bc 1300 // corr1->Write();
1301 // corr2->Write();
4e0af875 1302
5cd139bc 1303 for(Int_t iC=0;iC<nvarToPlot; iC++){
1304 for(Int_t iStep=0;iStep<3;iStep++){
1305 h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1306 }
1307 }
1308 file_projection->Close();
1309 for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1310 delete [] h;
6299d56c 1311
379592af 1312
379592af 1313}
1314
1315//___________________________________________________________________________
f2703bd2 1316void AliCFTaskVertexingHF::UserCreateOutputObjects()
1317{
5cd139bc 1318 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1319 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1320 //
1321 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
379592af 1322
5cd139bc 1323 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1324 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1325 AliPIDResponse *localPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
1326
1327 if (fCuts->GetIsUsePID() && fDecayChannel==22) {
1328
1329 fCuts->GetPidHF()->SetPidResponse(localPIDResponse);
5cd139bc 1330 fCuts->GetPidHF()->SetOldPid(kFALSE);
ea6f40e6 1331 AliRDHFCutsLctoV0* lcv0Cuts=dynamic_cast<AliRDHFCutsLctoV0*>(fCuts);
1332 if(lcv0Cuts){
1333 lcv0Cuts->GetPidV0pos()->SetPidResponse(localPIDResponse);
1334 lcv0Cuts->GetPidV0neg()->SetPidResponse(localPIDResponse);
1335 lcv0Cuts->GetPidV0pos()->SetOldPid(kFALSE);
1336 lcv0Cuts->GetPidV0neg()->SetOldPid(kFALSE);
1337 }
5cd139bc 1338 }
0e355b7f 1339
5cd139bc 1340 //slot #1
1341 OpenFile(1);
1342 const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
78471b33 1343 fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1344 fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1345 fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
5cd139bc 1346
1347 PostData(1,fHistEventsProcessed) ;
1348 PostData(2,fCFManager->GetParticleContainer()) ;
1349 PostData(3,fCorrelation) ;
0e355b7f 1350
379592af 1351}
1352
f2703bd2 1353
942a78ce 1354//_________________________________________________________________________
1355void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17a(){
1356 // ad-hoc weight function from ratio of
1357 // pt spectra from FONLL 2.76 TeV and
1358 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1359 if(fFuncWeight) delete fFuncWeight;
1360 fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1361 fFuncWeight->SetParameter(0,4.63891e-02);
1362 fFuncWeight->SetParameter(1,1.51674e+01);
1363 fFuncWeight->SetParameter(2,4.09941e-01);
96ca9111 1364 fUseWeight=kTRUE;
942a78ce 1365}
1366//_________________________________________________________________________
1367void AliCFTaskVertexingHF::SetPtWeightsFromDataPbPb276overLHC12a17a(){
1368 // ad-hoc weight function from ratio of
1369 // D0 pt spectra in PbPb 2011 0-10% centrality and
1370 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1371 if(fFuncWeight) delete fFuncWeight;
1372 fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1373 fFuncWeight->SetParameter(0,1.43116e-02);
1374 fFuncWeight->SetParameter(1,4.37758e+02);
1375 fFuncWeight->SetParameter(2,3.08583);
96ca9111 1376 fUseWeight=kTRUE;
942a78ce 1377}
4fa6d5f8 1378
1d9442de 1379//_________________________________________________________________________
1380void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17b(){
1381 // weight function from the ratio of the LHC12a17b MC
1382 // and FONLL calculations for pp data
1383 if(fFuncWeight) delete fFuncWeight;
1384 fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1385 fFuncWeight->SetParameters(1.92381e+01, 5.05055e+00, 1.05314e+01, 2.5, 1.88214e-03, 3.44871e+00, -9.74325e-02, 1.97671e+00, -3.21278e-01);
1386 fUseWeight=kTRUE;
1387}
1388
1389//_________________________________________________________________________
1390void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276andBAMPSoverLHC12a17b(){
1391 // weight function from the ratio of the LHC12a17b MC
1392 // and FONLL calculations for pp data
1393 // corrected by the BAMPS Raa calculation for 30-50% CC
1394 if(fFuncWeight) delete fFuncWeight;
1395 fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1396 fFuncWeight->SetParameters(6.10443e+00, 1.53487e+00, 1.99474e+00, 2.5, 5.51172e-03, 5.86590e+00, -5.46963e-01, 9.41201e-02, -1.64323e-01);
1397 fUseWeight=kTRUE;
1398}
1399
3825480a 1400//_________________________________________________________________________
1401void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC13d3(){
1402 // weight function from the ratio of the LHC12a17b MC
1403 // and FONLL calculations for pp data
1404 if(fFuncWeight) delete fFuncWeight;
1405 fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,30.);
1406 fFuncWeight->SetParameters(2.94999e+00,3.47032e+00,2.81278e+00,2.5,1.93370e-02,3.86865e+00,-1.54113e-01,8.86944e-02,2.56267e-02);
1407 fUseWeight=kTRUE;
1408}
1409
f2703bd2 1410//_________________________________________________________________________
1411Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1412{
5cd139bc 1413 //
1414 // calculating the weight to fill the container
1415 //
f2703bd2 1416
5cd139bc 1417 // FNOLL central:
1418 // p0 = 1.63297e-01 --> 0.322643
1419 // p1 = 2.96275e+00
1420 // p2 = 2.30301e+00
1421 // p3 = 2.50000e+00
1422
1423 // PYTHIA
1424 // p0 = 1.85906e-01 --> 0.36609
1425 // p1 = 1.94635e+00
1426 // p2 = 1.40463e+00
1427 // p3 = 2.50000e+00
1428
1429 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1430 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1431
1432 Double_t dndpt_func1 = dNdptFit(pt,func1);
1433 if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1434 Double_t dndpt_func2 = dNdptFit(pt,func2);
1435 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1436 return dndpt_func1/dndpt_func2;
f2703bd2 1437}
1438
1439//__________________________________________________________________________________________________
1440Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1441{
5cd139bc 1442 //
1443 // calculating dNdpt
1444 //
f2703bd2 1445
5cd139bc 1446 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1447 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
f2703bd2 1448
5cd139bc 1449 return dNdpt;
f2703bd2 1450}
6c62cb59 1451
17bf1a34 1452//__________________________________________________________________________________________________
1453Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1454 //
1455 // calculating the z-vtx weight for the given run range
1456 //
1457
1458 if(runnumber>146824 || runnumber<146803) return 1.0;
1459
1460 Double_t func1[3] = {1.0, -0.5, 6.5 };
1461 Double_t func2[3] = {1.0, -0.5, 5.5 };
1462
1463 Double_t dzFunc1 = DodzFit(z,func1);
1464 Double_t dzFunc2 = DodzFit(z,func2);
1465
1466 return dzFunc1/dzFunc2;
1467}
1468
1469//__________________________________________________________________________________________________
1470Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1471
1472 //
1473 // Gaussian z-vtx shape
1474 //
1475 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1476
1477 Double_t value = par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(z-par[1])*(z-par[1])/2./par[2]/par[2]);
1478
1479 return value;
1480}
6c62cb59 1481//__________________________________________________________________________________________________
afb24c40 1482Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1483 //
1484 // calculates the Nch weight using the measured and generateed Nch distributions
1485 //
1486 if(nch<=0) return 0.;
f26726ed 1487 if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
afb24c40 1488 Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1489 Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
4fa6d5f8 1490 Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1491 if(fUseTrackletsWeight) weight = pMC;
1492 return weight;
afb24c40 1493}
1494//__________________________________________________________________________________________________
1495void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1496 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
8e1c1209 1497 //
1498 // for Nch > 70 the points were obtained with a double NBD distribution fit
1499 // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
1500 // fit1->SetParameter(1,1.63); // k parameter
1501 // fit1->SetParameter(2,12.8); // mean multiplicity
1502 //
1503 Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
afb24c40 1504 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1505 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1506 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1507 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1508 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
8e1c1209 1509 60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1510 80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1511 100.50,102.50};
1512 Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
afb24c40 1513 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1514 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1515 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1516 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1517 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
8e1c1209 1518 0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1519 0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1520 0.00000258};
afb24c40 1521
1522 if(fHistoMeasNch) delete fHistoMeasNch;
8e1c1209 1523 fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1524 for(Int_t i=0; i<81; i++){
afb24c40 1525 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1526 fHistoMeasNch->SetBinError(i+1,0.);
1527 }
1528}
1529
5cd139bc 1530//__________________________________________________________________________________________________
6c62cb59 1531Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1532 // processes output of Ds is selected
1533 Bool_t keep=kFALSE;
1534 if(recoAnalysisCuts > 0){
1535 Int_t isKKpi=recoAnalysisCuts&1;
1536 Int_t ispiKK=recoAnalysisCuts&2;
1537 Int_t isPhiKKpi=recoAnalysisCuts&4;
1538 Int_t isPhipiKK=recoAnalysisCuts&8;
1539 Int_t isK0starKKpi=recoAnalysisCuts&16;
1540 Int_t isK0starpiKK=recoAnalysisCuts&32;
1541 if(fDsOption==1){
1542 if(isKKpi && isPhiKKpi) keep=kTRUE;
1543 if(ispiKK && isPhipiKK) keep=kTRUE;
1544 }
1545 else if(fDsOption==2){
1546 if(isKKpi && isK0starKKpi) keep=kTRUE;
1547 if(ispiKK && isK0starpiKK) keep=kTRUE;
1548 }
1549 else if(fDsOption==3)keep=kTRUE;
1550 }
1551 return keep;
1552}
5cd139bc 1553//__________________________________________________________________________________________________
1554Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1555 // processes output of Lc->V0+bnachelor is selected
1556 Bool_t keep=kFALSE;
1557 if(recoAnalysisCuts > 0){
1558
1559 Int_t isK0Sp = recoAnalysisCuts&1;
1560 Int_t isLambdaBarpi = recoAnalysisCuts&2;
1561 Int_t isLambdapi = recoAnalysisCuts&4;
1562
1563 if(fLctoV0bachelorOption==1){
1564 if(isK0Sp) keep=kTRUE;
1565 }
1566 else if(fLctoV0bachelorOption==2){
1567 if(isLambdaBarpi) keep=kTRUE;
1568 }
1569 else if(fLctoV0bachelorOption==4){
1570 if(isLambdapi) keep=kTRUE;
1571 }
1572 else if(fLctoV0bachelorOption==7) {
1573 if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1574 }
1575 }
1576 return keep;
1577}
f26726ed 1578
1579//____________________________________________________________________________
1580TProfile* AliCFTaskVertexingHF::GetEstimatorHistogram(const AliVEvent* event){
1581 // Get Estimator Histogram from period event->GetRunNumber();
1582 //
1583 // If you select SPD tracklets in |eta|<1 you should use type == 1
1584 //
1585
1586 Int_t runNo = event->GetRunNumber();
1587 Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1588 if(runNo>114930 && runNo<117223) period = 0;
1589 if(runNo>119158 && runNo<120830) period = 1;
1590 if(runNo>122373 && runNo<126438) period = 2;
1591 if(runNo>127711 && runNo<130841) period = 3;
1592 if(period<0 || period>3) return 0;
1593
1594 return fMultEstimatorAvg[period];
1595}