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