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