]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliCFTaskVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Reverted wrong commit (thanks Raoul for spotting it)
[u/mrichter/AliRoot.git] / PWG3 / 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>
38
39#include "AliCFTaskVertexingHF.h"
40#include "AliStack.h"
41#include "AliMCEvent.h"
42#include "AliCFManager.h"
43#include "AliCFContainer.h"
44#include "AliLog.h"
45#include "AliAnalysisManager.h"
46#include "AliAODHandler.h"
47#include "AliAODEvent.h"
48#include "AliAODRecoDecay.h"
49#include "AliAODRecoDecayHF.h"
50#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 51#include "AliAODRecoDecayHF3Prong.h"
52#include "AliAODRecoDecayHF4Prong.h"
53#include "AliAODRecoCascadeHF.h"
379592af 54#include "AliAODMCParticle.h"
55#include "AliAODMCHeader.h"
56#include "AliESDtrack.h"
57#include "TChain.h"
58#include "THnSparse.h"
59#include "TH2D.h"
60#include "AliESDtrackCuts.h"
f2703bd2 61#include "AliRDHFCuts.h"
379592af 62#include "AliRDHFCutsD0toKpi.h"
f2703bd2 63#include "AliRDHFCutsDplustoKpipi.h"
64#include "AliRDHFCutsDStartoKpipi.h"
65#include "AliRDHFCutsDstoKKpi.h"
66#include "AliRDHFCutsLctopKpi.h"
67#include "AliRDHFCutsD0toKpipipi.h"
379592af 68#include "AliCFVertexingHF2Prong.h"
4e600f58 69#include "AliCFVertexingHF3Prong.h"
379592af 70#include "AliCFVertexingHF.h"
f2703bd2 71#include "AliAnalysisDataSlot.h"
72#include "AliAnalysisDataContainer.h"
379592af 73
74//__________________________________________________________________________
75AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
f2703bd2 76 AliAnalysisTaskSE(),
f2703bd2 77 fCFManager(0x0),
78 fHistEventsProcessed(0x0),
79 fCorrelation(0x0),
80 fCountMC(0),
81 fCountAcc(0),
82 fCountVertex(0),
83 fCountRefit(0),
84 fCountReco(0),
85 fCountRecoAcc(0),
86 fCountRecoITSClusters(0),
87 fCountRecoPPR(0),
88 fCountRecoPID(0),
89 fEvents(0),
90 fDecayChannel(0),
91 fFillFromGenerated(kFALSE),
92 fOriginDselection(0),
93 fAcceptanceUnf(kTRUE),
94 fCuts(0),
95 fUseWeight(kFALSE),
96 fWeight(1.),
31da6b05 97 fNvar(0),
98 fPartName(""),
3ee5eb83 99 fDauNames(""),
100 fSign(2)
379592af 101{
102 //
103 //Default ctor
104 //
105}
106//___________________________________________________________________________
f2703bd2 107AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
108 AliAnalysisTaskSE(name),
f2703bd2 109 fCFManager(0x0),
110 fHistEventsProcessed(0x0),
111 fCorrelation(0x0),
112 fCountMC(0),
113 fCountAcc(0),
114 fCountVertex(0),
115 fCountRefit(0),
116 fCountReco(0),
117 fCountRecoAcc(0),
118 fCountRecoITSClusters(0),
119 fCountRecoPPR(0),
120 fCountRecoPID(0),
121 fEvents(0),
122 fDecayChannel(0),
123 fFillFromGenerated(kFALSE),
124 fOriginDselection(0),
125 fAcceptanceUnf(kTRUE),
126 fCuts(cuts),
127 fUseWeight(kFALSE),
128 fWeight(1.),
31da6b05 129 fNvar(0),
130 fPartName(""),
3ee5eb83 131 fDauNames(""),
132 fSign(2)
379592af 133{
134 //
135 // Constructor. Initialization of Inputs and Outputs
136 //
379592af 137 /*
f2703bd2 138 DefineInput(0) and DefineOutput(0)
139 are taken care of by AliAnalysisTaskSE constructor
140 */
379592af 141 DefineOutput(1,TH1I::Class());
142 DefineOutput(2,AliCFContainer::Class());
143 DefineOutput(3,THnSparseD::Class());
f2703bd2 144 DefineOutput(4,AliRDHFCuts::Class());
379592af 145
146 fCuts->PrintAll();
147}
148
149//___________________________________________________________________________
150AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
151{
152 //
153 // Assignment operator
154 //
155 if (this!=&c) {
156 AliAnalysisTaskSE::operator=(c) ;
379592af 157 fCFManager = c.fCFManager;
158 fHistEventsProcessed = c.fHistEventsProcessed;
159 fCuts = c.fCuts;
160 }
161 return *this;
162}
163
164//___________________________________________________________________________
165AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 166 AliAnalysisTaskSE(c),
f2703bd2 167 fCFManager(c.fCFManager),
168 fHistEventsProcessed(c.fHistEventsProcessed),
169 fCorrelation(c.fCorrelation),
170 fCountMC(c.fCountMC),
171 fCountAcc(c.fCountAcc),
172 fCountVertex(c.fCountVertex),
173 fCountRefit(c.fCountRefit),
174 fCountReco(c.fCountReco),
175 fCountRecoAcc(c.fCountRecoAcc),
176 fCountRecoITSClusters(c.fCountRecoITSClusters),
177 fCountRecoPPR(c.fCountRecoPPR),
178 fCountRecoPID(c.fCountRecoPID),
179 fEvents(c.fEvents),
180 fDecayChannel(c.fDecayChannel),
181 fFillFromGenerated(c.fFillFromGenerated),
182 fOriginDselection(c.fOriginDselection),
183 fAcceptanceUnf(c.fAcceptanceUnf),
184 fCuts(c.fCuts),
185 fUseWeight(c.fUseWeight),
186 fWeight(c.fWeight),
31da6b05 187 fNvar(c.fNvar),
188 fPartName(c.fPartName),
3ee5eb83 189 fDauNames(c.fDauNames),
190 fSign(c.fSign)
379592af 191{
192 //
193 // Copy Constructor
194 //
195}
196
197//___________________________________________________________________________
f2703bd2 198AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
199{
379592af 200 //
201 //destructor
202 //
203 if (fCFManager) delete fCFManager ;
204 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 205 if (fCorrelation) delete fCorrelation ;
379592af 206 if (fCuts) delete fCuts;
207}
208
f2703bd2 209//_________________________________________________________________________-
210void AliCFTaskVertexingHF::Init()
211{
212 //
213 // Initialization
214 //
215
216 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
217 AliRDHFCuts *copyfCuts = 0x0;
8a983038 218 if (!fCuts){
219 AliFatal("No cuts defined - Exiting...");
220 return;
221 }
222
f2703bd2 223 switch (fDecayChannel){
224 case 2:{
225 copyfCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fCuts)));
226 fNvar = 13;
31da6b05 227 fPartName="D0";
228 fDauNames="K+pi";
f2703bd2 229 break;
230 }
231 case 21:{
232 copyfCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
233 fNvar = 13;
31da6b05 234 fPartName="Dstar";
235 fDauNames="K+pi+pi";
f2703bd2 236 break;
237 }
238 case 31:{
239 copyfCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
240 fNvar = 12;
31da6b05 241 fPartName="Dplus";
242 fDauNames="K+pi+pi";
f2703bd2 243 break;
244 }
245 case 32:{
246 copyfCuts = new AliRDHFCutsLctopKpi(*(dynamic_cast<AliRDHFCutsLctopKpi*>(fCuts)));
31da6b05 247 fNvar = 12;
248 fPartName="Lambdac";
249 fDauNames="p+K+pi";
f2703bd2 250 break;
251 }
252 case 33:{
253 copyfCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
31da6b05 254 fNvar = 12;
255 fPartName="Ds";
256 fDauNames="K+K+pi";
f2703bd2 257 break;
258 }
259 case 4:{
260 copyfCuts = new AliRDHFCutsD0toKpipipi(*(dynamic_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
261 fNvar = 13;
31da6b05 262 fPartName="D0";
263 fDauNames="K+pi+pi+pi";
f2703bd2 264 break;
265 }
266 default:
267 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
268 break;
269 }
270
271 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
8a983038 272 if (copyfCuts){
273 copyfCuts->SetName(nameoutput);
274
275 //Post the data
276 PostData(4, copyfCuts);
277 }
278 else{
279 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
280 }
f2703bd2 281
282 return;
283}
284
379592af 285//_________________________________________________
286void AliCFTaskVertexingHF::UserExec(Option_t *)
287{
379592af 288 //
289 // Main loop function
290 //
291
292 PostData(1,fHistEventsProcessed) ;
293 PostData(2,fCFManager->GetParticleContainer()) ;
294 PostData(3,fCorrelation) ;
295
296 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts();
297
298 if (fFillFromGenerated){
299 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
300 }
301
302 if (!fInputEvent) {
303 Error("UserExec","NO EVENT FOUND!");
304 return;
305 }
306
307 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
308
f2703bd2 309 TClonesArray *arrayBranch=0;
379592af 310
311 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 312 // In case there is an AOD handler writing a standard AOD, use the AOD
313 // event in memory rather than the input (ESD) event.
314 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
315 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
316 // have to taken from the AOD event hold by the AliAODExtension
317 AliAODHandler* aodHandler = (AliAODHandler*)
318 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
319 if(aodHandler->GetExtensions()) {
320 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
321 AliAODEvent *aodFromExt = ext->GetAOD();
322
323 switch (fDecayChannel){
324 case 2:{
325 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
326 break;
327 }
328 case 21:{
329 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
330 break;
331 }
332 case 31:
333 case 32:
334 case 33:{
335 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
336 break;
337 }
338 case 4:{
339 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
340 break;
341 }
342 default:
343 break;
344 }
345 }
346 }
347 else {
348 switch (fDecayChannel){
349 case 2:{
350 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
351 break;
352 }
353 case 21:{
354 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
355 break;
356 }
357 case 31:
358 case 32:
359 case 33:{
360 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
361 break;
362 }
363 case 4:{
364 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
365 break;
366 }
367 default:
368 break;
369 }
379592af 370 }
371
372 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
373 if (!aodVtx) return;
f2703bd2 374
375 if (!arrayBranch) {
376 AliError("Could not find array of HF vertices");
377 return;
379592af 378 }
379
380 fEvents++;
6e2e4f50 381
379592af 382 fCFManager->SetRecEventInfo(aodEvent);
383 fCFManager->SetMCEventInfo(aodEvent);
384
379592af 385 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 386
379592af 387 //loop on the MC event
388
389 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
390 if (!mcArray) {
391 AliError("Could not find Monte-Carlo in AOD");
392 return;
393 }
394 Int_t icountMC = 0;
395 Int_t icountAcc = 0;
396 Int_t icountReco = 0;
397 Int_t icountVertex = 0;
398 Int_t icountRefit = 0;
399 Int_t icountRecoAcc = 0;
400 Int_t icountRecoITSClusters = 0;
401 Int_t icountRecoPPR = 0;
402 Int_t icountRecoPID = 0;
379592af 403 Int_t cquarks = 0;
f2703bd2 404
379592af 405 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
406 if (!mcHeader) {
407 AliError("Could not find MC Header in AOD");
408 return;
409 }
e11ae259 410
411 Double_t* containerInput = new Double_t[fNvar];
412 Double_t* containerInputMC = new Double_t[fNvar];
413
f2703bd2 414
415 AliCFVertexingHF* cfVtxHF=0x0;
416 switch (fDecayChannel){
417 case 2:{
6e2e4f50 418 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
f2703bd2 419 break;
420 }
421 case 21:{
422 // cfVtxHF = new AliCFVertexingHFCascade(mcArray, originDselection); // not there yet
423 break;
424 }
425 case 31:
426 case 32:
427 case 33:{
4e600f58 428 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
f2703bd2 429 break;
430 }
431 case 4:{
432 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
433 break;
434 }
435 default:
436 break;
437 }
8a983038 438 if (!cfVtxHF){
439 AliError("No AliCFVertexingHF initialized");
440 return;
441 }
f2703bd2 442
379592af 443 Double_t zPrimVertex = aodVtx ->GetZ();
444 Double_t zMCVertex = mcHeader->GetVtxZ();
f2703bd2 445
379592af 446 //General settings: vertex, feed down and fill reco container with generated values.
447 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
448 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
449 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
f2703bd2 450
379592af 451 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
379592af 452
f2703bd2 453 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
8a983038 454 if (!mcPart){
455 AliError("Failed casting particle from MC array!, Skipping particle");
456 continue;
457 }
f2703bd2 458 // check the MC-level cuts, must be the desidered particle
459 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
460
461 cfVtxHF->SetMCCandidateParam(iPart);
462 cfVtxHF->SetNVar(fNvar);
463 //counting c quarks
464 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
465
6e2e4f50 466 if (!(cfVtxHF->SetLabelArray())){
467 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
468 continue;
469 }
470
f2703bd2 471 //check the candiate family at MC level
472 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
6e2e4f50 473 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 474 continue;
475 }
476 else{
6e2e4f50 477 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 478 }
479
480 //Fill the MC container
481 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 482 if (mcContainerFilled) {
483 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
484 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
485 //MC Limited Acceptance
486 if (TMath::Abs(containerInputMC[1]) < 0.5) {
487 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
488 AliDebug(3,"MC Lim Acc container filled\n");
489 }
490
491 //MC
492 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
493 icountMC++;
494 AliDebug(3,"MC cointainer filled \n");
495
496 // MC in acceptance
497 // check the MC-Acceptance level cuts
498 // since standard CF functions are not applicable, using Kine Cuts on daughters
499 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
500 if (mcAccepStep){
501 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
502 AliDebug(3,"MC acceptance cut passed\n");
503 icountAcc++;
504
505 //MC Vertex step
506 if (fCuts->IsEventSelected(aodEvent)){
507 // filling the container if the vertex is ok
508 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
509 AliDebug(3,"Vertex cut passed and container filled\n");
510 icountVertex++;
511
512 //mc Refit requirement
513 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
514 if (mcRefitStep){
515 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
516 AliDebug(3,"MC Refit cut passed and container filled\n");
517 icountRefit++;
518 }
519 else{
520 AliDebug(3,"MC Refit cut not passed\n");
521 continue;
522 }
523 }
524 else{
379592af 525 AliDebug (3, "MC vertex step not passed\n");
526 continue;
f2703bd2 527 }
528 }
529 else{
530 AliDebug (3, "MC in acceptance step not passed\n");
531 continue;
532 }
533 }
534 else {
535 AliDebug (3, "MC container not filled\n");
379592af 536 }
537 }
f2703bd2 538
539 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 540 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
541 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
542 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
543 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
f2703bd2 544
379592af 545 // Now go to rec level
546 fCountMC += icountMC;
547 fCountAcc += icountAcc;
548 fCountVertex+= icountVertex;
549 fCountRefit+= icountRefit;
f2703bd2 550
6e2e4f50 551 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 552
f2703bd2 553 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 554 AliAODRecoDecayHF* charmCandidate=0x0;
555 switch (fDecayChannel){
556 case 2:{
557 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
558 break;
559 }
560 case 21:{
561 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
562 break;
563 }
564 case 31:
565 case 32:
566 case 33:{
567 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
568 break;
569 }
570 case 4:{
571 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
572 break;
573 }
574 default:
575 break;
576 }
577
379592af 578 Bool_t unsetvtx=kFALSE;
f2703bd2 579 if(!charmCandidate->GetOwnPrimaryVtx()) {
580 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 581 unsetvtx=kTRUE;
582 }
f2703bd2 583
584 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 585 if (!signAssociation){
f2703bd2 586 charmCandidate = 0x0;
379592af 587 continue;
588 }
3ee5eb83 589
590 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
591 if (isPartOrAntipart == 0){
592 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
593 continue;
594 }
595
379592af 596 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
597 if (recoContFilled){
f2703bd2 598
379592af 599 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 600
379592af 601 //Reco Step
602 Bool_t recoStep = cfVtxHF->RecoStep();
603 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 604
379592af 605 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 606 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 607 icountReco++;
f2703bd2 608 AliDebug(3,"Reco step passed and container filled\n");
609
379592af 610 //Reco in the acceptance -- take care of UNFOLDING!!!!
611 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
612 if (recoAcceptanceStep) {
f2703bd2 613 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 614 icountRecoAcc++;
f2703bd2 615 AliDebug(3,"Reco acceptance cut passed and container filled\n");
616
379592af 617 if(fAcceptanceUnf){
618 Double_t fill[4]; //fill response matrix
619 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
620 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 621 }
622
379592af 623 //Number of ITS cluster requirements
f2703bd2 624 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 625 if (recoITSnCluster){
f2703bd2 626 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 627 icountRecoITSClusters++;
f2703bd2 628 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
629
630 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
631 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
632 fCuts->SetUsePID(kFALSE);
633 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
634
635 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
636 if (recoAnalysisCuts >= 8){
637 recoAnalysisCuts -= 8; // removing K0star mass
638 }
639 if (recoAnalysisCuts >= 4){
640 recoAnalysisCuts -= 4; // removing Phi mass
641 }
642 }
379592af 643
f2703bd2 644 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
645 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
646 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
647 icountRecoPPR++;
648 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
649 //pid selection
650 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
651 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
652 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
653 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
654 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
655 icountRecoPID++;
656 AliDebug(3,"Reco PID cuts passed and container filled \n");
657 if(!fAcceptanceUnf){
658 Double_t fill[4]; //fill response matrix
659 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
660 if (bUnfolding) fCorrelation->Fill(fill);
661 }
662 }
663 else {
664 AliDebug(3, "Analysis Cuts step not passed \n");
665 continue;
666 }
379592af 667 }
668 else {
f2703bd2 669 AliDebug(3, "PID selection not passed \n");
670 continue;
379592af 671 }
672 }
673 else{
f2703bd2 674 AliDebug(3, "Number of ITS cluster step not passed\n");
675 continue;
379592af 676 }
677 }
678 else{
f2703bd2 679 AliDebug(3, "Reco acceptance step not passed\n");
680 continue;
379592af 681 }
682 }
683 else {
f2703bd2 684 AliDebug(3, "Reco step not passed\n");
685 continue;
379592af 686 }
687 }
688
f2703bd2 689 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
690 } // end loop on candidate
691
379592af 692 fCountReco+= icountReco;
693 fCountRecoAcc+= icountRecoAcc;
694 fCountRecoITSClusters+= icountRecoITSClusters;
695 fCountRecoPPR+= icountRecoPPR;
696 fCountRecoPID+= icountRecoPID;
697
698 fHistEventsProcessed->Fill(0);
f2703bd2 699
700 delete[] containerInput;
701 containerInput = 0x0;
702 delete[] containerInputMC;
703 containerInputMC = 0x0;
704 delete cfVtxHF;
705
379592af 706}
707
708//___________________________________________________________________________
709void AliCFTaskVertexingHF::Terminate(Option_t*)
710{
711 // The Terminate() function is the last function to be called during
712 // a query. It always runs on the client, it can be used to present
713 // the results graphically or save the results to file.
714
715 AliAnalysisTaskSE::Terminate();
31da6b05 716
717 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
718 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
719 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));
720 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));
721 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
722 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));
723 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));
724 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));
725 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 726
727 // draw some example plots....
728
729 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
730 if(!cont) {
731 printf("CONTAINER NOT FOUND\n");
732 return;
733 }
734 // projecting the containers to obtain histograms
735 // first argument = variable, second argument = step
736
31da6b05 737 TH1D* h[3][12];
738 for(Int_t iC=0;iC<12; iC++){
739 // MC-level
740 h[0][iC] = cont->ShowProjection(iC,0);
741 // MC-Acceptance level
742 h[1][iC] = cont->ShowProjection(iC,1);
743 // Reco-level
744 h[2][iC] = cont->ShowProjection(iC,4);
745 }
746 TString titles[12];
747 if(fDecayChannel==31){
748 titles[0]="pT_Dplus (GeV/c)";
749 titles[1]="rapidity";
750 titles[2]="phi (rad)";
751 titles[3]="cT (#mum)";
752 titles[4]="cosPointingAngle";
753 titles[5]="pT_1 (GeV/c)";
754 titles[6]="pT_2 (GeV/c)";
755 titles[7]="pT_3 (GeV/c)";
756 titles[8]="d0_1 (#mum)";
757 titles[9]="d0_2 (#mum)";
758 titles[10]="d0_3 (#mum)";
759 titles[11]="zVertex (cm)";
760 }else{
761 titles[0]="pT_D0 (GeV/c)";
762 titles[1]="rapidity";
763 titles[2]="cosThetaStar";
764 titles[3]="pT_pi (GeV/c)";
765 titles[4]="pT_K (Gev/c)";
766 titles[5]="cT (#mum)";
767 titles[6]="dca (#mum)";
768 titles[7]="d0_pi (#mum)";
769 titles[8]="d0_K (#mum)";
770 titles[9]="d0xd0 (#mum^2)";
771 titles[10]="cosPointingAngle";
772 titles[11]="phi (rad)";
773 }
774 Int_t markers[12]={20,24,21,25,27,28,
775 20,24,21,25,27,28};
776 Int_t colors[3]={2,8,4};
777 for(Int_t iC=0;iC<12; iC++){
778 for(Int_t iStep=0;iStep<3;iStep++){
779 h[iStep][iC]->SetTitle(titles[iC].Data());
780 h[iStep][iC]->GetXaxis()->SetTitle(titles[iC].Data());
781 Double_t maxh=h[iStep][iC]->GetMaximum();
782 h[iStep][iC]->GetYaxis()->SetRangeUser(0,maxh*1.2);
783 h[iStep][iC]->SetMarkerStyle(markers[iC]);
784 h[iStep][iC]->SetMarkerColor(colors[iStep]);
785 }
786 }
787
379592af 788
379592af 789
790 gStyle->SetCanvasColor(0);
791 gStyle->SetFrameFillColor(0);
792 gStyle->SetTitleFillColor(0);
793 gStyle->SetStatColor(0);
794
795 // drawing in 2 separate canvas for a matter of clearity
31da6b05 796 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
379592af 797 c1->Divide(3,3);
31da6b05 798 Int_t iPad=1;
799 for(Int_t iVar=0; iVar<3; iVar++){
800 c1->cd(iPad++);
801 h[0][iVar]->Draw("p");
802 c1->cd(iPad++);
803 h[1][iVar]->Draw("p");
804 c1->cd(iPad++);
805 h[2][iVar]->Draw("p");
806 }
379592af 807
31da6b05 808 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
379592af 809 c2->Divide(3,3);
31da6b05 810 iPad=1;
811 for(Int_t iVar=3; iVar<6; iVar++){
812 c2->cd(iPad++);
813 h[0][iVar]->Draw("p");
814 c2->cd(iPad++);
815 h[1][iVar]->Draw("p");
816 c2->cd(iPad++);
817 h[2][iVar]->Draw("p");
818 }
819
379592af 820
31da6b05 821 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
379592af 822 c3->Divide(3,3);
31da6b05 823 iPad=1;
824 for(Int_t iVar=6; iVar<9; iVar++){
825 c3->cd(iPad++);
826 h[0][iVar]->Draw("p");
827 c3->cd(iPad++);
828 h[1][iVar]->Draw("p");
829 c3->cd(iPad++);
830 h[2][iVar]->Draw("p");
831 }
832
379592af 833
31da6b05 834 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
379592af 835 c4->Divide(3,3);
31da6b05 836 iPad=1;
837 for(Int_t iVar=9; iVar<11; iVar++){
838 c4->cd(iPad++);
839 h[0][iVar]->Draw("p");
840 c4->cd(iPad++);
841 h[1][iVar]->Draw("p");
842 c4->cd(iPad++);
843 h[2][iVar]->Draw("p");
844 }
845
379592af 846
847 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
848
849 TH2D* corr1 =hcorr->Projection(0,2);
850 TH2D* corr2 = hcorr->Projection(1,3);
851
852 TCanvas * c7 =new TCanvas("c7New","",800,400);
853 c7->Divide(2,1);
854 c7->cd(1);
855 corr1->Draw("text");
856 c7->cd(2);
857 corr2->Draw("text");
858
859
860 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
861
862 corr1->Write();
863 corr2->Write();
31da6b05 864 for(Int_t iC=0;iC<12; iC++){
865 for(Int_t iStep=0;iStep<3;iStep++){
866 h[iStep][iC]->Write(Form("Step%d_%s",iStep,titles[iC].Data()));
867 }
868 }
379592af 869 file_projection->Close();
870
871
872
873 /*
874 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
875 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
876 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
877 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
878
879 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
880 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
881 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
882 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
883 */
884}
885
886//___________________________________________________________________________
f2703bd2 887void AliCFTaskVertexingHF::UserCreateOutputObjects()
888{
379592af 889 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
890 //TO BE SET BEFORE THE EXECUTION OF THE TASK
891 //
892 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
893
894 //slot #1
895 OpenFile(1);
896 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
897}
898
f2703bd2 899
900//_________________________________________________________________________
901Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
902{
903 //
904 // calculating the weight to fill the container
905 //
906
907 // FNOLL central:
908 // p0 = 1.63297e-01 --> 0.322643
909 // p1 = 2.96275e+00
910 // p2 = 2.30301e+00
911 // p3 = 2.50000e+00
912
913 // PYTHIA
914 // p0 = 1.85906e-01 --> 0.36609
915 // p1 = 1.94635e+00
916 // p2 = 1.40463e+00
917 // p3 = 2.50000e+00
918
919 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
920 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
921
922 Double_t dndpt_func1 = dNdptFit(pt,func1);
923 Double_t dndpt_func2 = dNdptFit(pt,func2);
924 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
925 return dndpt_func1/dndpt_func2;
926}
927
928//__________________________________________________________________________________________________
929Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
930{
931 //
932 // calculating dNdpt
933 //
934
935 Double_t denom = TMath::Power((pt/par[1]), par[3] );
936 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
937
938 return dNdpt;
939}