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