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