Updates in PbPb cuts (Andrea Rossi)
[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"
367e9aa3 70#include "AliCFVertexingHFCascade.h"
379592af 71#include "AliCFVertexingHF.h"
f2703bd2 72#include "AliAnalysisDataSlot.h"
73#include "AliAnalysisDataContainer.h"
379592af 74
75//__________________________________________________________________________
76AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
f2703bd2 77 AliAnalysisTaskSE(),
f2703bd2 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.),
31da6b05 98 fNvar(0),
99 fPartName(""),
3ee5eb83 100 fDauNames(""),
b7af2639 101 fSign(2),
fbec9fa9 102 fCentralitySelection(kTRUE),
0ada222f 103 fFakeSelection(0),
104 fRejectIfNoQuark(kTRUE),
6c62cb59 105 fUseMCVertex(kFALSE),
106 fDsOption(1)
379592af 107{
108 //
109 //Default ctor
110 //
111}
112//___________________________________________________________________________
f2703bd2 113AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
114 AliAnalysisTaskSE(name),
f2703bd2 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.),
31da6b05 135 fNvar(0),
136 fPartName(""),
3ee5eb83 137 fDauNames(""),
b7af2639 138 fSign(2),
fbec9fa9 139 fCentralitySelection(kTRUE),
0ada222f 140 fFakeSelection(0),
141 fRejectIfNoQuark(kTRUE),
6c62cb59 142 fUseMCVertex(kFALSE),
143 fDsOption(1)
379592af 144{
145 //
146 // Constructor. Initialization of Inputs and Outputs
147 //
379592af 148 /*
f2703bd2 149 DefineInput(0) and DefineOutput(0)
150 are taken care of by AliAnalysisTaskSE constructor
151 */
379592af 152 DefineOutput(1,TH1I::Class());
153 DefineOutput(2,AliCFContainer::Class());
154 DefineOutput(3,THnSparseD::Class());
f2703bd2 155 DefineOutput(4,AliRDHFCuts::Class());
379592af 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) ;
379592af 168 fCFManager = c.fCFManager;
169 fHistEventsProcessed = c.fHistEventsProcessed;
170 fCuts = c.fCuts;
171 }
172 return *this;
173}
174
175//___________________________________________________________________________
176AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 177 AliAnalysisTaskSE(c),
f2703bd2 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),
31da6b05 198 fNvar(c.fNvar),
199 fPartName(c.fPartName),
3ee5eb83 200 fDauNames(c.fDauNames),
b7af2639 201 fSign(c.fSign),
fbec9fa9 202 fCentralitySelection(c.fCentralitySelection),
0ada222f 203 fFakeSelection(c.fFakeSelection),
204 fRejectIfNoQuark(c.fRejectIfNoQuark),
6c62cb59 205 fUseMCVertex(c.fUseMCVertex),
206 fDsOption(c.fDsOption)
379592af 207{
208 //
209 // Copy Constructor
210 //
211}
212
213//___________________________________________________________________________
f2703bd2 214AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
215{
379592af 216 //
217 //destructor
218 //
219 if (fCFManager) delete fCFManager ;
220 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 221 if (fCorrelation) delete fCorrelation ;
379592af 222 if (fCuts) delete fCuts;
223}
224
f2703bd2 225//_________________________________________________________________________-
226void AliCFTaskVertexingHF::Init()
227{
228 //
229 // Initialization
230 //
231
232 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
233 AliRDHFCuts *copyfCuts = 0x0;
8a983038 234 if (!fCuts){
235 AliFatal("No cuts defined - Exiting...");
236 return;
237 }
238
f2703bd2 239 switch (fDecayChannel){
240 case 2:{
314ec0fe 241 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
1f780958 242 fNvar = 15;
31da6b05 243 fPartName="D0";
244 fDauNames="K+pi";
f2703bd2 245 break;
246 }
247 case 21:{
314ec0fe 248 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
1f780958 249 fNvar = 15;
31da6b05 250 fPartName="Dstar";
251 fDauNames="K+pi+pi";
f2703bd2 252 break;
253 }
254 case 31:{
314ec0fe 255 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
1f780958 256 fNvar = 14;
31da6b05 257 fPartName="Dplus";
258 fDauNames="K+pi+pi";
f2703bd2 259 break;
260 }
261 case 32:{
314ec0fe 262 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
1c44b663 263 fNvar = 18;
31da6b05 264 fPartName="Lambdac";
265 fDauNames="p+K+pi";
f2703bd2 266 break;
267 }
268 case 33:{
314ec0fe 269 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
1f780958 270 fNvar = 14;
31da6b05 271 fPartName="Ds";
272 fDauNames="K+K+pi";
f2703bd2 273 break;
274 }
275 case 4:{
314ec0fe 276 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
1f780958 277 fNvar = 15;
31da6b05 278 fPartName="D0";
279 fDauNames="K+pi+pi+pi";
f2703bd2 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();
8a983038 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 }
f2703bd2 297
298 return;
299}
300
379592af 301//_________________________________________________
302void AliCFTaskVertexingHF::UserExec(Option_t *)
303{
379592af 304 //
305 // Main loop function
306 //
307
308 PostData(1,fHistEventsProcessed) ;
309 PostData(2,fCFManager->GetParticleContainer()) ;
0e355b7f 310 PostData(3,fCorrelation) ;
b7af2639 311
379592af 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
f2703bd2 323 TClonesArray *arrayBranch=0;
379592af 324
325 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 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 }
379592af 384 }
385
386 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
387 if (!aodVtx) return;
f2703bd2 388
389 if (!arrayBranch) {
390 AliError("Could not find array of HF vertices");
391 return;
379592af 392 }
393
394 fEvents++;
6e2e4f50 395
379592af 396 fCFManager->SetRecEventInfo(aodEvent);
397 fCFManager->SetMCEventInfo(aodEvent);
398
379592af 399 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 400
379592af 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;
379592af 417 Int_t cquarks = 0;
f2703bd2 418
379592af 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 }
e11ae259 424
425 Double_t* containerInput = new Double_t[fNvar];
426 Double_t* containerInputMC = new Double_t[fNvar];
427
f2703bd2 428
429 AliCFVertexingHF* cfVtxHF=0x0;
430 switch (fDecayChannel){
431 case 2:{
367e9aa3 432 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
433 break;
f2703bd2 434 }
435 case 21:{
367e9aa3 436 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
437 break;
f2703bd2 438 }
439 case 31:
440 case 32:
441 case 33:{
4e600f58 442 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
f2703bd2 443 break;
444 }
445 case 4:{
446 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
447 break;
448 }
449 default:
450 break;
451 }
8a983038 452 if (!cfVtxHF){
453 AliError("No AliCFVertexingHF initialized");
d49092e3 454 delete[] containerInput;
455 delete[] containerInputMC;
8a983038 456 return;
457 }
f2703bd2 458
379592af 459 Double_t zPrimVertex = aodVtx ->GetZ();
460 Double_t zMCVertex = mcHeader->GetVtxZ();
751b8274 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()));
3655ecf2 463 delete[] containerInput;
464 delete[] containerInputMC;
751b8274 465 return;
466 }
467
2bf2e62b 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 }
367e9aa3 474 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
2bf2e62b 475 }
476 else {
477 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
478 trackCuts[iProng]=fCuts->GetTrackCuts();
479 }
480 }
481
379592af 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);
b7af2639 486 cfVtxHF->SetNVar(fNvar);
fbec9fa9 487 cfVtxHF->SetFakeSelection(fFakeSelection);
0ada222f 488 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
b7af2639 489
af548a24 490 // switch-off the trigger class selection (doesn't work for MC)
491 fCuts->SetTriggerClass("");
492
0ada222f 493 // MC vertex, to be used, in case, for pp
494 if (fUseMCVertex) fCuts->SetUseMCVertex();
af548a24 495
72156e2e 496 if (fCentralitySelection){ // keep only the requested centrality
661d6c29 497 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
498 delete[] containerInput;
499 delete[] containerInputMC;
500 delete [] trackCuts;
501 return;
502 }
72156e2e 503 } else { // keep all centralities
504 fCuts->SetMinCentrality(0.);
505 fCuts->SetMaxCentrality(100.);
8b1b055d 506 }
507
b7af2639 508
509 Float_t centValue = fCuts->GetCentrality(aodEvent);
510 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 511
379592af 512 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
b7af2639 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 }
2bf2e62b 531
b7af2639 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 }
f2703bd2 537 else{
6e2e4f50 538 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 539 }
540
541 //Fill the MC container
542 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 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
2bf2e62b 574 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
f2703bd2 575 if (mcRefitStep){
576 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
577 AliDebug(3,"MC Refit cut passed and container filled\n");
578 icountRefit++;
b7af2639 579
580
f2703bd2 581 }
582 else{
583 AliDebug(3,"MC Refit cut not passed\n");
584 continue;
585 }
586 }
587 else{
379592af 588 AliDebug (3, "MC vertex step not passed\n");
589 continue;
f2703bd2 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");
379592af 599 }
600 }
f2703bd2 601
602 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 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()));
f2703bd2 607
379592af 608 // Now go to rec level
609 fCountMC += icountMC;
610 fCountAcc += icountAcc;
611 fCountVertex+= icountVertex;
612 fCountRefit+= icountRefit;
f2703bd2 613
6e2e4f50 614 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 615
f2703bd2 616 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 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
379592af 641 Bool_t unsetvtx=kFALSE;
f2703bd2 642 if(!charmCandidate->GetOwnPrimaryVtx()) {
643 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 644 unsetvtx=kTRUE;
645 }
f2703bd2 646
647 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 648 if (!signAssociation){
f2703bd2 649 charmCandidate = 0x0;
379592af 650 continue;
651 }
3ee5eb83 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
6694a2a8 659
379592af 660 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
661 if (recoContFilled){
6694a2a8 662
663 // weight according to pt
664 if (fUseWeight)fWeight = GetWeight(containerInput[0]);
665
379592af 666 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 667
379592af 668 //Reco Step
669 Bool_t recoStep = cfVtxHF->RecoStep();
670 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 671
379592af 672 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 673 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 674 icountReco++;
f2703bd2 675 AliDebug(3,"Reco step passed and container filled\n");
676
379592af 677 //Reco in the acceptance -- take care of UNFOLDING!!!!
2bf2e62b 678 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
379592af 679 if (recoAcceptanceStep) {
f2703bd2 680 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 681 icountRecoAcc++;
f2703bd2 682 AliDebug(3,"Reco acceptance cut passed and container filled\n");
683
379592af 684 if(fAcceptanceUnf){
685 Double_t fill[4]; //fill response matrix
686 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
687 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 688 }
689
379592af 690 //Number of ITS cluster requirements
f2703bd2 691 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 692 if (recoITSnCluster){
f2703bd2 693 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 694 icountRecoITSClusters++;
f2703bd2 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
6c62cb59 702 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
703 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
704 if(keepDs) recoAnalysisCuts=3;
f2703bd2 705 }
6c62cb59 706
379592af 707
f2703bd2 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);
6c62cb59 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
f2703bd2 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 }
379592af 737 }
738 else {
f2703bd2 739 AliDebug(3, "PID selection not passed \n");
740 continue;
379592af 741 }
742 }
743 else{
f2703bd2 744 AliDebug(3, "Number of ITS cluster step not passed\n");
745 continue;
379592af 746 }
747 }
748 else{
f2703bd2 749 AliDebug(3, "Reco acceptance step not passed\n");
750 continue;
379592af 751 }
752 }
753 else {
f2703bd2 754 AliDebug(3, "Reco step not passed\n");
755 continue;
379592af 756 }
757 }
758
f2703bd2 759 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
760 } // end loop on candidate
761
379592af 762 fCountReco+= icountReco;
763 fCountRecoAcc+= icountRecoAcc;
764 fCountRecoITSClusters+= icountRecoITSClusters;
765 fCountRecoPPR+= icountRecoPPR;
766 fCountRecoPID+= icountRecoPID;
767
768 fHistEventsProcessed->Fill(0);
f2703bd2 769
770 delete[] containerInput;
f2703bd2 771 delete[] containerInputMC;
f2703bd2 772 delete cfVtxHF;
2bf2e62b 773 if (trackCuts){
774 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
775 // delete [] trackCuts[i];
776 // }
777 delete [] trackCuts;
778 }
379592af 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();
31da6b05 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));
379592af 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
31da6b05 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
379592af 861
379592af 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
31da6b05 869 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
379592af 870 c1->Divide(3,3);
31da6b05 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 }
379592af 880
31da6b05 881 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
379592af 882 c2->Divide(3,3);
31da6b05 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
379592af 893
31da6b05 894 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
379592af 895 c3->Divide(3,3);
31da6b05 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
379592af 906
31da6b05 907 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
379592af 908 c4->Divide(3,3);
31da6b05 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
379592af 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();
31da6b05 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 }
379592af 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//___________________________________________________________________________
f2703bd2 960void AliCFTaskVertexingHF::UserCreateOutputObjects()
961{
379592af 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) ;
0e355b7f 970
971 PostData(1,fHistEventsProcessed) ;
972 PostData(2,fCFManager->GetParticleContainer()) ;
973 PostData(3,fCorrelation) ;
974
379592af 975}
976
f2703bd2 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}
6c62cb59 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}