]> 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
#94948: Change in the TOF response and porting into the v5-02 release
[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>
35#include <TH1I.h>
36#include <TStyle.h>
37#include <TFile.h>
5b37c6e5 38#include <TF1.h>
379592af 39
40#include "AliCFTaskVertexingHF.h"
41#include "AliStack.h"
42#include "AliMCEvent.h"
43#include "AliCFManager.h"
44#include "AliCFContainer.h"
45#include "AliLog.h"
46#include "AliAnalysisManager.h"
47#include "AliAODHandler.h"
48#include "AliAODEvent.h"
49#include "AliAODRecoDecay.h"
50#include "AliAODRecoDecayHF.h"
51#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 52#include "AliAODRecoDecayHF3Prong.h"
53#include "AliAODRecoDecayHF4Prong.h"
54#include "AliAODRecoCascadeHF.h"
379592af 55#include "AliAODMCParticle.h"
56#include "AliAODMCHeader.h"
57#include "AliESDtrack.h"
58#include "TChain.h"
59#include "THnSparse.h"
60#include "TH2D.h"
61#include "AliESDtrackCuts.h"
f2703bd2 62#include "AliRDHFCuts.h"
379592af 63#include "AliRDHFCutsD0toKpi.h"
f2703bd2 64#include "AliRDHFCutsDplustoKpipi.h"
65#include "AliRDHFCutsDStartoKpipi.h"
66#include "AliRDHFCutsDstoKKpi.h"
67#include "AliRDHFCutsLctopKpi.h"
68#include "AliRDHFCutsD0toKpipipi.h"
379592af 69#include "AliCFVertexingHF2Prong.h"
4e600f58 70#include "AliCFVertexingHF3Prong.h"
367e9aa3 71#include "AliCFVertexingHFCascade.h"
379592af 72#include "AliCFVertexingHF.h"
f2703bd2 73#include "AliAnalysisDataSlot.h"
74#include "AliAnalysisDataContainer.h"
379592af 75
76//__________________________________________________________________________
77AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
f2703bd2 78 AliAnalysisTaskSE(),
f2703bd2 79 fCFManager(0x0),
80 fHistEventsProcessed(0x0),
81 fCorrelation(0x0),
82 fCountMC(0),
83 fCountAcc(0),
84 fCountVertex(0),
85 fCountRefit(0),
86 fCountReco(0),
87 fCountRecoAcc(0),
88 fCountRecoITSClusters(0),
89 fCountRecoPPR(0),
90 fCountRecoPID(0),
91 fEvents(0),
92 fDecayChannel(0),
93 fFillFromGenerated(kFALSE),
94 fOriginDselection(0),
95 fAcceptanceUnf(kTRUE),
96 fCuts(0),
97 fUseWeight(kFALSE),
98 fWeight(1.),
17bf1a34 99 fUseFlatPtWeight(kFALSE),
100 fUseZWeight(kFALSE),
31da6b05 101 fNvar(0),
102 fPartName(""),
3ee5eb83 103 fDauNames(""),
b7af2639 104 fSign(2),
fbec9fa9 105 fCentralitySelection(kTRUE),
0ada222f 106 fFakeSelection(0),
107 fRejectIfNoQuark(kTRUE),
6c62cb59 108 fUseMCVertex(kFALSE),
ec5288c3 109 fDsOption(1),
c83eda41 110 fGenDsOption(3),
5b37c6e5 111 fConfiguration(kCheetah), // by default, setting the fast configuration
112 fFuncWeight(0x0)
379592af 113{
114 //
115 //Default ctor
116 //
117}
118//___________________________________________________________________________
5b37c6e5 119AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
f2703bd2 120 AliAnalysisTaskSE(name),
f2703bd2 121 fCFManager(0x0),
122 fHistEventsProcessed(0x0),
123 fCorrelation(0x0),
124 fCountMC(0),
125 fCountAcc(0),
126 fCountVertex(0),
127 fCountRefit(0),
128 fCountReco(0),
129 fCountRecoAcc(0),
130 fCountRecoITSClusters(0),
131 fCountRecoPPR(0),
132 fCountRecoPID(0),
133 fEvents(0),
134 fDecayChannel(0),
135 fFillFromGenerated(kFALSE),
136 fOriginDselection(0),
137 fAcceptanceUnf(kTRUE),
138 fCuts(cuts),
139 fUseWeight(kFALSE),
140 fWeight(1.),
17bf1a34 141 fUseFlatPtWeight(kFALSE),
142 fUseZWeight(kFALSE),
31da6b05 143 fNvar(0),
144 fPartName(""),
3ee5eb83 145 fDauNames(""),
b7af2639 146 fSign(2),
fbec9fa9 147 fCentralitySelection(kTRUE),
0ada222f 148 fFakeSelection(0),
149 fRejectIfNoQuark(kTRUE),
6c62cb59 150 fUseMCVertex(kFALSE),
ec5288c3 151 fDsOption(1),
c83eda41 152 fGenDsOption(3),
5b37c6e5 153 fConfiguration(kCheetah), // by default, setting the fast configuration
154 fFuncWeight(func)
379592af 155{
156 //
157 // Constructor. Initialization of Inputs and Outputs
158 //
379592af 159 /*
f2703bd2 160 DefineInput(0) and DefineOutput(0)
161 are taken care of by AliAnalysisTaskSE constructor
162 */
379592af 163 DefineOutput(1,TH1I::Class());
164 DefineOutput(2,AliCFContainer::Class());
165 DefineOutput(3,THnSparseD::Class());
f2703bd2 166 DefineOutput(4,AliRDHFCuts::Class());
379592af 167
168 fCuts->PrintAll();
169}
170
171//___________________________________________________________________________
172AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
173{
174 //
175 // Assignment operator
176 //
177 if (this!=&c) {
178 AliAnalysisTaskSE::operator=(c) ;
379592af 179 fCFManager = c.fCFManager;
180 fHistEventsProcessed = c.fHistEventsProcessed;
181 fCuts = c.fCuts;
5b37c6e5 182 fFuncWeight = c.fFuncWeight;
379592af 183 }
184 return *this;
185}
186
187//___________________________________________________________________________
188AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 189 AliAnalysisTaskSE(c),
f2703bd2 190 fCFManager(c.fCFManager),
191 fHistEventsProcessed(c.fHistEventsProcessed),
192 fCorrelation(c.fCorrelation),
193 fCountMC(c.fCountMC),
194 fCountAcc(c.fCountAcc),
195 fCountVertex(c.fCountVertex),
196 fCountRefit(c.fCountRefit),
197 fCountReco(c.fCountReco),
198 fCountRecoAcc(c.fCountRecoAcc),
199 fCountRecoITSClusters(c.fCountRecoITSClusters),
200 fCountRecoPPR(c.fCountRecoPPR),
201 fCountRecoPID(c.fCountRecoPID),
202 fEvents(c.fEvents),
203 fDecayChannel(c.fDecayChannel),
204 fFillFromGenerated(c.fFillFromGenerated),
205 fOriginDselection(c.fOriginDselection),
206 fAcceptanceUnf(c.fAcceptanceUnf),
207 fCuts(c.fCuts),
208 fUseWeight(c.fUseWeight),
209 fWeight(c.fWeight),
17bf1a34 210 fUseFlatPtWeight(c.fUseFlatPtWeight),
211 fUseZWeight(c.fUseZWeight),
31da6b05 212 fNvar(c.fNvar),
213 fPartName(c.fPartName),
3ee5eb83 214 fDauNames(c.fDauNames),
b7af2639 215 fSign(c.fSign),
fbec9fa9 216 fCentralitySelection(c.fCentralitySelection),
0ada222f 217 fFakeSelection(c.fFakeSelection),
218 fRejectIfNoQuark(c.fRejectIfNoQuark),
6c62cb59 219 fUseMCVertex(c.fUseMCVertex),
ec5288c3 220 fDsOption(c.fDsOption),
c83eda41 221 fGenDsOption(c.fGenDsOption),
5b37c6e5 222 fConfiguration(c.fConfiguration),
223 fFuncWeight(c.fFuncWeight)
379592af 224{
225 //
226 // Copy Constructor
227 //
228}
229
230//___________________________________________________________________________
f2703bd2 231AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
232{
379592af 233 //
234 //destructor
235 //
236 if (fCFManager) delete fCFManager ;
237 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 238 if (fCorrelation) delete fCorrelation ;
379592af 239 if (fCuts) delete fCuts;
5b37c6e5 240 if (fFuncWeight) delete fFuncWeight;
379592af 241}
242
f2703bd2 243//_________________________________________________________________________-
244void AliCFTaskVertexingHF::Init()
245{
246 //
247 // Initialization
248 //
249
250 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
17bf1a34 251 if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
f2703bd2 252 AliRDHFCuts *copyfCuts = 0x0;
8a983038 253 if (!fCuts){
254 AliFatal("No cuts defined - Exiting...");
255 return;
256 }
257
f2703bd2 258 switch (fDecayChannel){
259 case 2:{
314ec0fe 260 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
ec5288c3 261 switch (fConfiguration) {
262 case kSnail: // slow configuration: all variables in
263 fNvar = 16;
264 break;
265 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
266 fNvar = 8;
267 break;
268 }
31da6b05 269 fPartName="D0";
270 fDauNames="K+pi";
f2703bd2 271 break;
272 }
273 case 21:{
314ec0fe 274 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
ec5288c3 275 switch (fConfiguration) {
276 case kSnail: // slow configuration: all variables in
277 fNvar = 16;
278 break;
279 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
280 fNvar = 8;
281 break;
282 }
31da6b05 283 fPartName="Dstar";
284 fDauNames="K+pi+pi";
f2703bd2 285 break;
286 }
287 case 31:{
314ec0fe 288 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
ec5288c3 289 switch (fConfiguration) {
290 case kSnail: // slow configuration: all variables in
291 fNvar = 14;
292 break;
293 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
294 fNvar = 8;
295 break;
296 }
31da6b05 297 fPartName="Dplus";
298 fDauNames="K+pi+pi";
f2703bd2 299 break;
300 }
301 case 32:{
314ec0fe 302 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
ec5288c3 303 switch (fConfiguration) {
304 case kSnail: // slow configuration: all variables in
305 fNvar = 18;
306 break;
307 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
308 fNvar = 8;
309 break;
310 }
31da6b05 311 fPartName="Lambdac";
312 fDauNames="p+K+pi";
f2703bd2 313 break;
314 }
315 case 33:{
314ec0fe 316 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
ec5288c3 317 switch (fConfiguration) {
318 case kSnail: // slow configuration: all variables in
319 fNvar = 14;
320 break;
321 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
322 fNvar = 8;
323 break;
324 }
31da6b05 325 fPartName="Ds";
326 fDauNames="K+K+pi";
f2703bd2 327 break;
328 }
329 case 4:{
314ec0fe 330 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
ec5288c3 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 }
31da6b05 339 fPartName="D0";
340 fDauNames="K+pi+pi+pi";
f2703bd2 341 break;
342 }
343 default:
344 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
345 break;
346 }
347
348 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
8a983038 349 if (copyfCuts){
350 copyfCuts->SetName(nameoutput);
351
352 //Post the data
353 PostData(4, copyfCuts);
354 }
355 else{
356 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
357 }
f2703bd2 358
359 return;
360}
361
379592af 362//_________________________________________________
363void AliCFTaskVertexingHF::UserExec(Option_t *)
364{
379592af 365 //
366 // Main loop function
367 //
368
369 PostData(1,fHistEventsProcessed) ;
370 PostData(2,fCFManager->GetParticleContainer()) ;
0e355b7f 371 PostData(3,fCorrelation) ;
b7af2639 372
379592af 373 if (fFillFromGenerated){
374 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
375 }
376
377 if (!fInputEvent) {
378 Error("UserExec","NO EVENT FOUND!");
379 return;
380 }
381
382 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
383
f2703bd2 384 TClonesArray *arrayBranch=0;
379592af 385
386 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 387 // In case there is an AOD handler writing a standard AOD, use the AOD
388 // event in memory rather than the input (ESD) event.
389 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
390 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
391 // have to taken from the AOD event hold by the AliAODExtension
392 AliAODHandler* aodHandler = (AliAODHandler*)
393 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
394 if(aodHandler->GetExtensions()) {
395 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
396 AliAODEvent *aodFromExt = ext->GetAOD();
397
398 switch (fDecayChannel){
399 case 2:{
400 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
401 break;
402 }
403 case 21:{
404 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
405 break;
406 }
407 case 31:
408 case 32:
409 case 33:{
410 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
411 break;
412 }
413 case 4:{
414 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
415 break;
416 }
417 default:
418 break;
419 }
420 }
421 }
422 else {
423 switch (fDecayChannel){
424 case 2:{
425 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
426 break;
427 }
428 case 21:{
429 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
430 break;
431 }
432 case 31:
433 case 32:
434 case 33:{
435 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
436 break;
437 }
438 case 4:{
439 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
440 break;
441 }
442 default:
443 break;
444 }
379592af 445 }
ec5288c3 446
379592af 447 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
448 if (!aodVtx) return;
f2703bd2 449
450 if (!arrayBranch) {
451 AliError("Could not find array of HF vertices");
452 return;
379592af 453 }
454
455 fEvents++;
6e2e4f50 456
379592af 457 fCFManager->SetRecEventInfo(aodEvent);
458 fCFManager->SetMCEventInfo(aodEvent);
459
379592af 460 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 461
379592af 462 //loop on the MC event
463
464 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
465 if (!mcArray) {
466 AliError("Could not find Monte-Carlo in AOD");
467 return;
468 }
469 Int_t icountMC = 0;
470 Int_t icountAcc = 0;
471 Int_t icountReco = 0;
472 Int_t icountVertex = 0;
473 Int_t icountRefit = 0;
474 Int_t icountRecoAcc = 0;
475 Int_t icountRecoITSClusters = 0;
476 Int_t icountRecoPPR = 0;
477 Int_t icountRecoPID = 0;
379592af 478 Int_t cquarks = 0;
f2703bd2 479
379592af 480 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
481 if (!mcHeader) {
482 AliError("Could not find MC Header in AOD");
483 return;
484 }
e11ae259 485
486 Double_t* containerInput = new Double_t[fNvar];
487 Double_t* containerInputMC = new Double_t[fNvar];
488
f2703bd2 489
490 AliCFVertexingHF* cfVtxHF=0x0;
491 switch (fDecayChannel){
492 case 2:{
367e9aa3 493 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
494 break;
f2703bd2 495 }
496 case 21:{
367e9aa3 497 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
498 break;
f2703bd2 499 }
500 case 31:
501 case 32:
502 case 33:{
4e600f58 503 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
c83eda41 504 if(fDecayChannel==33){
505 cfVtxHF->SetGeneratedDsOption(fGenDsOption);
506 }
f2703bd2 507 break;
508 }
509 case 4:{
510 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
511 break;
512 }
513 default:
514 break;
515 }
8a983038 516 if (!cfVtxHF){
517 AliError("No AliCFVertexingHF initialized");
d49092e3 518 delete[] containerInput;
519 delete[] containerInputMC;
8a983038 520 return;
521 }
f2703bd2 522
379592af 523 Double_t zPrimVertex = aodVtx ->GetZ();
524 Double_t zMCVertex = mcHeader->GetVtxZ();
17bf1a34 525 Int_t runnumber = aodEvent->GetRunNumber();
526 if(fUseZWeight) fWeight = GetZWeight(zMCVertex,runnumber);
527
751b8274 528 if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
529 AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
3655ecf2 530 delete[] containerInput;
531 delete[] containerInputMC;
751b8274 532 return;
533 }
534
2bf2e62b 535 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
536 if (fDecayChannel == 21){
537 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
538 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
539 trackCuts[iProng]=fCuts->GetTrackCuts();
540 }
367e9aa3 541 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
2bf2e62b 542 }
543 else {
544 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
545 trackCuts[iProng]=fCuts->GetTrackCuts();
546 }
547 }
548
379592af 549 //General settings: vertex, feed down and fill reco container with generated values.
550 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
551 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
552 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
b7af2639 553 cfVtxHF->SetNVar(fNvar);
fbec9fa9 554 cfVtxHF->SetFakeSelection(fFakeSelection);
0ada222f 555 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
ec5288c3 556 cfVtxHF->SetConfiguration(fConfiguration);
b7af2639 557
af548a24 558 // switch-off the trigger class selection (doesn't work for MC)
559 fCuts->SetTriggerClass("");
560
0ada222f 561 // MC vertex, to be used, in case, for pp
562 if (fUseMCVertex) fCuts->SetUseMCVertex();
af548a24 563
72156e2e 564 if (fCentralitySelection){ // keep only the requested centrality
661d6c29 565 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
566 delete[] containerInput;
567 delete[] containerInputMC;
568 delete [] trackCuts;
569 return;
570 }
72156e2e 571 } else { // keep all centralities
572 fCuts->SetMinCentrality(0.);
573 fCuts->SetMaxCentrality(100.);
8b1b055d 574 }
575
b7af2639 576 Float_t centValue = fCuts->GetCentrality(aodEvent);
577 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 578
ec5288c3 579 // number of tracklets - multiplicity estimator
580 Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles
581 cfVtxHF->SetMultiplicity(multiplicity);
582
379592af 583 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
b7af2639 584 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
585 if (!mcPart){
586 AliError("Failed casting particle from MC array!, Skipping particle");
587 continue;
588 }
589 // check the MC-level cuts, must be the desidered particle
590 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
591 continue; // 0 stands for MC level
592 }
593 cfVtxHF->SetMCCandidateParam(iPart);
594
595 //counting c quarks
596 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
597
598 if (!(cfVtxHF->SetLabelArray())){
599 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
600 continue;
601 }
2bf2e62b 602
b7af2639 603 //check the candiate family at MC level
604 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
605 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
606 continue;
607 }
f2703bd2 608 else{
6e2e4f50 609 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 610 }
611
612 //Fill the MC container
613 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 614 if (mcContainerFilled) {
5b37c6e5 615 if (fUseWeight){
616 if (fFuncWeight){ // using user-defined function
617 AliDebug(2,"Using function");
618 fWeight = fFuncWeight->Eval(containerInputMC[0]);
619 }
620 else{ // using FONLL
621 AliDebug(2,"Using FONLL");
622 fWeight = GetWeight(containerInputMC[0]);
623 }
624 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
625 }
f2703bd2 626 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
627 //MC Limited Acceptance
628 if (TMath::Abs(containerInputMC[1]) < 0.5) {
629 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
630 AliDebug(3,"MC Lim Acc container filled\n");
631 }
632
633 //MC
634 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
635 icountMC++;
636 AliDebug(3,"MC cointainer filled \n");
637
638 // MC in acceptance
639 // check the MC-Acceptance level cuts
640 // since standard CF functions are not applicable, using Kine Cuts on daughters
641 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
642 if (mcAccepStep){
643 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
644 AliDebug(3,"MC acceptance cut passed\n");
645 icountAcc++;
646
647 //MC Vertex step
648 if (fCuts->IsEventSelected(aodEvent)){
649 // filling the container if the vertex is ok
650 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
651 AliDebug(3,"Vertex cut passed and container filled\n");
652 icountVertex++;
653
654 //mc Refit requirement
2bf2e62b 655 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
f2703bd2 656 if (mcRefitStep){
657 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
658 AliDebug(3,"MC Refit cut passed and container filled\n");
659 icountRefit++;
660 }
661 else{
662 AliDebug(3,"MC Refit cut not passed\n");
663 continue;
664 }
665 }
666 else{
379592af 667 AliDebug (3, "MC vertex step not passed\n");
668 continue;
f2703bd2 669 }
670 }
671 else{
672 AliDebug (3, "MC in acceptance step not passed\n");
673 continue;
674 }
675 }
676 else {
677 AliDebug (3, "MC container not filled\n");
379592af 678 }
679 }
f2703bd2 680
681 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 682 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
683 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
684 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
685 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
f2703bd2 686
379592af 687 // Now go to rec level
688 fCountMC += icountMC;
689 fCountAcc += icountAcc;
690 fCountVertex+= icountVertex;
691 fCountRefit+= icountRefit;
f2703bd2 692
6e2e4f50 693 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 694
f2703bd2 695 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 696 AliAODRecoDecayHF* charmCandidate=0x0;
697 switch (fDecayChannel){
698 case 2:{
699 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
700 break;
701 }
702 case 21:{
703 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
704 break;
705 }
706 case 31:
707 case 32:
708 case 33:{
709 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
710 break;
711 }
712 case 4:{
713 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
714 break;
715 }
716 default:
717 break;
718 }
719
379592af 720 Bool_t unsetvtx=kFALSE;
f2703bd2 721 if(!charmCandidate->GetOwnPrimaryVtx()) {
722 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 723 unsetvtx=kTRUE;
724 }
f2703bd2 725
726 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 727 if (!signAssociation){
f2703bd2 728 charmCandidate = 0x0;
379592af 729 continue;
730 }
3ee5eb83 731
732 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
733 if (isPartOrAntipart == 0){
734 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
735 continue;
736 }
737
6694a2a8 738
379592af 739 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
740 if (recoContFilled){
6694a2a8 741
742 // weight according to pt
5b37c6e5 743 if (fUseWeight){
744 if (fFuncWeight){ // using user-defined function
745 AliDebug(2, "Using function");
746 fWeight = fFuncWeight->Eval(containerInput[0]);
747 }
748 else{ // using FONLL
749 AliDebug(2, "Using FONLL");
750 fWeight = GetWeight(containerInput[0]);
751 }
752 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
753 }
6694a2a8 754
379592af 755 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 756
379592af 757 //Reco Step
758 Bool_t recoStep = cfVtxHF->RecoStep();
759 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 760
379592af 761 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 762 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 763 icountReco++;
f2703bd2 764 AliDebug(3,"Reco step passed and container filled\n");
765
379592af 766 //Reco in the acceptance -- take care of UNFOLDING!!!!
2bf2e62b 767 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
379592af 768 if (recoAcceptanceStep) {
f2703bd2 769 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 770 icountRecoAcc++;
f2703bd2 771 AliDebug(3,"Reco acceptance cut passed and container filled\n");
772
379592af 773 if(fAcceptanceUnf){
774 Double_t fill[4]; //fill response matrix
775 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
776 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 777 }
778
379592af 779 //Number of ITS cluster requirements
f2703bd2 780 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 781 if (recoITSnCluster){
f2703bd2 782 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 783 icountRecoITSClusters++;
f2703bd2 784 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
785
786 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
787 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
788 fCuts->SetUsePID(kFALSE);
789 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
790
6c62cb59 791 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
792 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
793 if(keepDs) recoAnalysisCuts=3;
f2703bd2 794 }
6c62cb59 795
379592af 796
f2703bd2 797 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
798 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
799 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
800 icountRecoPPR++;
801 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
802 //pid selection
803 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
804 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
805 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
6c62cb59 806
807 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
808 Bool_t keepDs=ProcessDs(recoPidSelection);
809 if(keepDs) recoPidSelection=3;
810 }
811
f2703bd2 812 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
813 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
814 icountRecoPID++;
815 AliDebug(3,"Reco PID cuts passed and container filled \n");
816 if(!fAcceptanceUnf){
817 Double_t fill[4]; //fill response matrix
818 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
819 if (bUnfolding) fCorrelation->Fill(fill);
820 }
821 }
822 else {
823 AliDebug(3, "Analysis Cuts step not passed \n");
824 continue;
825 }
379592af 826 }
827 else {
f2703bd2 828 AliDebug(3, "PID selection not passed \n");
829 continue;
379592af 830 }
831 }
832 else{
f2703bd2 833 AliDebug(3, "Number of ITS cluster step not passed\n");
834 continue;
379592af 835 }
836 }
837 else{
f2703bd2 838 AliDebug(3, "Reco acceptance step not passed\n");
839 continue;
379592af 840 }
841 }
842 else {
f2703bd2 843 AliDebug(3, "Reco step not passed\n");
844 continue;
379592af 845 }
846 }
847
f2703bd2 848 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
849 } // end loop on candidate
850
379592af 851 fCountReco+= icountReco;
852 fCountRecoAcc+= icountRecoAcc;
853 fCountRecoITSClusters+= icountRecoITSClusters;
854 fCountRecoPPR+= icountRecoPPR;
855 fCountRecoPID+= icountRecoPID;
856
857 fHistEventsProcessed->Fill(0);
f2703bd2 858
859 delete[] containerInput;
f2703bd2 860 delete[] containerInputMC;
f2703bd2 861 delete cfVtxHF;
2bf2e62b 862 if (trackCuts){
863 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
864 // delete [] trackCuts[i];
865 // }
866 delete [] trackCuts;
867 }
ec5288c3 868
869
379592af 870}
871
872//___________________________________________________________________________
873void AliCFTaskVertexingHF::Terminate(Option_t*)
874{
875 // The Terminate() function is the last function to be called during
876 // a query. It always runs on the client, it can be used to present
877 // the results graphically or save the results to file.
878
879 AliAnalysisTaskSE::Terminate();
31da6b05 880
881 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
882 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
883 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));
884 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));
885 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
886 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));
887 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));
888 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));
889 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 890
891 // draw some example plots....
ec5288c3 892 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
379592af 893 if(!cont) {
894 printf("CONTAINER NOT FOUND\n");
895 return;
896 }
897 // projecting the containers to obtain histograms
898 // first argument = variable, second argument = step
ec5288c3 899
900 TH1D** h = new TH1D*[3];
901 if (fConfiguration == kSnail){
902 //h = new TH1D[3][12];
903 for (Int_t ih = 0; ih<3; ih++){
904 h[ih] = new TH1D[12];
905 }
906 for(Int_t iC=1;iC<12; iC++){
907 // MC-level
908 h[0][iC] = *(cont->ShowProjection(iC,0));
909 // MC-Acceptance level
910 h[1][iC] = *(cont->ShowProjection(iC,1));
911 // Reco-level
912 h[2][iC] = *(cont->ShowProjection(iC,4));
913 }
31da6b05 914 }
ec5288c3 915 else{
916 //h = new TH1D[3][12];
917 for (Int_t ih = 0; ih<3; ih++){
918 h[ih] = new TH1D[8];
919 }
920 for(Int_t iC=0;iC<8; iC++){
921 // MC-level
922 h[0][iC] = *(cont->ShowProjection(iC,0));
923 // MC-Acceptance level
924 h[1][iC] = *(cont->ShowProjection(iC,1));
925 // Reco-level
926 h[2][iC] = *(cont->ShowProjection(iC,4));
927 }
928 }
ec5288c3 929 TString* titles;
930 Int_t nvarToPlot = 0;
931 if (fConfiguration == kSnail){
932 nvarToPlot = 12;
933 titles = new TString[nvarToPlot];
934 if(fDecayChannel==31){
935 titles[0]="pT_Dplus (GeV/c)";
936 titles[1]="rapidity";
937 titles[2]="phi (rad)";
938 titles[3]="cT (#mum)";
939 titles[4]="cosPointingAngle";
940 titles[5]="pT_1 (GeV/c)";
941 titles[6]="pT_2 (GeV/c)";
942 titles[7]="pT_3 (GeV/c)";
943 titles[8]="d0_1 (#mum)";
944 titles[9]="d0_2 (#mum)";
945 titles[10]="d0_3 (#mum)";
946 titles[11]="zVertex (cm)";
947 }else{
948 titles[0]="pT_D0 (GeV/c)";
949 titles[1]="rapidity";
950 titles[2]="cosThetaStar";
951 titles[3]="pT_pi (GeV/c)";
952 titles[4]="pT_K (Gev/c)";
953 titles[5]="cT (#mum)";
954 titles[6]="dca (#mum)";
955 titles[7]="d0_pi (#mum)";
956 titles[8]="d0_K (#mum)";
957 titles[9]="d0xd0 (#mum^2)";
958 titles[10]="cosPointingAngle";
959 titles[11]="phi (rad)";
960
961 }
962 }
963 else{
964 nvarToPlot = 8;
965 titles = new TString[nvarToPlot];
966 titles[0]="pT_candidate (GeV/c)";
967 titles[1]="rapidity";
968 titles[2]="cT (#mum)";
969 titles[3]="phi";
970 titles[4]="z_{vtx}";
971 titles[5]="centrality";
972 titles[6]="fake";
973 titles[7]="multiplicity";
31da6b05 974 }
ec5288c3 975
31da6b05 976 Int_t markers[12]={20,24,21,25,27,28,
977 20,24,21,25,27,28};
978 Int_t colors[3]={2,8,4};
ec5288c3 979 for(Int_t iC=0;iC<nvarToPlot; iC++){
980 for(Int_t iStep=0;iStep<3;iStep++){
981 h[iStep][iC].SetTitle(titles[iC].Data());
982 h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
983 Double_t maxh=h[iStep][iC].GetMaximum();
984 h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
985 h[iStep][iC].SetMarkerStyle(markers[iC]);
986 h[iStep][iC].SetMarkerColor(colors[iStep]);
987 }
31da6b05 988 }
379592af 989
990 gStyle->SetCanvasColor(0);
991 gStyle->SetFrameFillColor(0);
992 gStyle->SetTitleFillColor(0);
993 gStyle->SetStatColor(0);
994
995 // drawing in 2 separate canvas for a matter of clearity
ec5288c3 996 TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
997 c1->Divide(3,4);
31da6b05 998 Int_t iPad=1;
ec5288c3 999 for(Int_t iVar=0; iVar<4; iVar++){
1000 c1->cd(iPad++);
5b37c6e5 1001 h[0][iVar].DrawCopy("p");
ec5288c3 1002 c1->cd(iPad++);
5b37c6e5 1003 h[1][iVar].DrawCopy("p");
ec5288c3 1004 c1->cd(iPad++);
5b37c6e5 1005 h[2][iVar].DrawCopy("p");
31da6b05 1006 }
379592af 1007
ec5288c3 1008 TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1009 c2->Divide(3,4);
31da6b05 1010 iPad=1;
ec5288c3 1011 for(Int_t iVar=4; iVar<8; iVar++){
1012 c2->cd(iPad++);
5b37c6e5 1013 h[0][iVar].DrawCopy("p");
ec5288c3 1014 c2->cd(iPad++);
5b37c6e5 1015 h[1][iVar].DrawCopy("p");
ec5288c3 1016 c2->cd(iPad++);
5b37c6e5 1017 h[2][iVar].DrawCopy("p");
31da6b05 1018 }
1019
ec5288c3 1020 if (fConfiguration == kSnail){
1021 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1022 c3->Divide(3,4);
1023 iPad=1;
1024 for(Int_t iVar=8; iVar<12; iVar++){
1025 c3->cd(iPad++);
5b37c6e5 1026 h[0][iVar].DrawCopy("p");
ec5288c3 1027 c3->cd(iPad++);
5b37c6e5 1028 h[1][iVar].DrawCopy("p");
ec5288c3 1029 c3->cd(iPad++);
5b37c6e5 1030 h[2][iVar].DrawCopy("p");
ec5288c3 1031 }
31da6b05 1032 }
1033
379592af 1034
1035 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1036
1037 TH2D* corr1 =hcorr->Projection(0,2);
1038 TH2D* corr2 = hcorr->Projection(1,3);
1039
ec5288c3 1040 TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
379592af 1041 c7->Divide(2,1);
1042 c7->cd(1);
5b37c6e5 1043 corr1->DrawCopy("text");
379592af 1044 c7->cd(2);
5b37c6e5 1045 corr2->DrawCopy("text");
379592af 1046
379592af 1047 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1048
1049 corr1->Write();
1050 corr2->Write();
ec5288c3 1051 for(Int_t iC=0;iC<nvarToPlot; iC++){
1052 for(Int_t iStep=0;iStep<3;iStep++){
1053 h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1054 }
31da6b05 1055 }
379592af 1056 file_projection->Close();
5b37c6e5 1057 for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
dc818294 1058 delete [] h;
6299d56c 1059
379592af 1060
379592af 1061}
1062
1063//___________________________________________________________________________
f2703bd2 1064void AliCFTaskVertexingHF::UserCreateOutputObjects()
1065{
379592af 1066 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1067 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1068 //
1069 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1070
1071 //slot #1
1072 OpenFile(1);
1073 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
0e355b7f 1074
1075 PostData(1,fHistEventsProcessed) ;
1076 PostData(2,fCFManager->GetParticleContainer()) ;
1077 PostData(3,fCorrelation) ;
1078
379592af 1079}
1080
f2703bd2 1081
1082//_________________________________________________________________________
1083Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1084{
1085 //
1086 // calculating the weight to fill the container
1087 //
1088
1089 // FNOLL central:
1090 // p0 = 1.63297e-01 --> 0.322643
1091 // p1 = 2.96275e+00
1092 // p2 = 2.30301e+00
1093 // p3 = 2.50000e+00
1094
1095 // PYTHIA
1096 // p0 = 1.85906e-01 --> 0.36609
1097 // p1 = 1.94635e+00
1098 // p2 = 1.40463e+00
1099 // p3 = 2.50000e+00
1100
1101 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1102 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1103
1104 Double_t dndpt_func1 = dNdptFit(pt,func1);
17bf1a34 1105 if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
f2703bd2 1106 Double_t dndpt_func2 = dNdptFit(pt,func2);
1107 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1108 return dndpt_func1/dndpt_func2;
1109}
1110
1111//__________________________________________________________________________________________________
1112Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1113{
1114 //
1115 // calculating dNdpt
1116 //
1117
1118 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1119 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1120
1121 return dNdpt;
1122}
6c62cb59 1123
17bf1a34 1124//__________________________________________________________________________________________________
1125Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1126 //
1127 // calculating the z-vtx weight for the given run range
1128 //
1129
1130 if(runnumber>146824 || runnumber<146803) return 1.0;
1131
1132 Double_t func1[3] = {1.0, -0.5, 6.5 };
1133 Double_t func2[3] = {1.0, -0.5, 5.5 };
1134
1135 Double_t dzFunc1 = DodzFit(z,func1);
1136 Double_t dzFunc2 = DodzFit(z,func2);
1137
1138 return dzFunc1/dzFunc2;
1139}
1140
1141//__________________________________________________________________________________________________
1142Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1143
1144 //
1145 // Gaussian z-vtx shape
1146 //
1147 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1148
1149 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]);
1150
1151 return value;
1152}
6c62cb59 1153
1154//__________________________________________________________________________________________________
1155Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1156 // processes output of Ds is selected
1157 Bool_t keep=kFALSE;
1158 if(recoAnalysisCuts > 0){
1159 Int_t isKKpi=recoAnalysisCuts&1;
1160 Int_t ispiKK=recoAnalysisCuts&2;
1161 Int_t isPhiKKpi=recoAnalysisCuts&4;
1162 Int_t isPhipiKK=recoAnalysisCuts&8;
1163 Int_t isK0starKKpi=recoAnalysisCuts&16;
1164 Int_t isK0starpiKK=recoAnalysisCuts&32;
1165 if(fDsOption==1){
1166 if(isKKpi && isPhiKKpi) keep=kTRUE;
1167 if(ispiKK && isPhipiKK) keep=kTRUE;
1168 }
1169 else if(fDsOption==2){
1170 if(isKKpi && isK0starKKpi) keep=kTRUE;
1171 if(ispiKK && isK0starpiKK) keep=kTRUE;
1172 }
1173 else if(fDsOption==3)keep=kTRUE;
1174 }
1175 return keep;
1176}