]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
allow for selection of signal range for high gain channels (useful e.g. for cosmic...
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
CommitLineData
f5ec6c30 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// 11 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25//-----------------------------------------------------------------------
a5d92cb4 26//-----------------------------------------------------------------------
27// Base class for HF Unfolding (pt and eta)
28// correlation matrix filled at Acceptance and PPR level
29// Author: A.Grelli , Utrecht - agrelli@uu.nl
30//-----------------------------------------------------------------------
f5ec6c30 31#include <TCanvas.h>
32#include <TParticle.h>
c180f65d 33#include <TDatabasePDG.h>
f5ec6c30 34#include <TH1I.h>
35#include <TStyle.h>
844d235c 36#include <TFile.h>
f5ec6c30 37
38#include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
39#include "AliStack.h"
40#include "AliMCEvent.h"
41#include "AliCFManager.h"
42#include "AliCFContainer.h"
43#include "AliLog.h"
44#include "AliAODEvent.h"
45#include "AliAODRecoDecay.h"
46#include "AliAODRecoDecayHF.h"
47#include "AliAODRecoDecayHF2Prong.h"
48#include "AliAODMCParticle.h"
aed46c56 49#include "AliESDtrack.h"
a5d92cb4 50#include "TChain.h"
51#include "THnSparse.h"
52#include "TH2D.h"
f5ec6c30 53//__________________________________________________________________________
54AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
55 AliAnalysisTaskSE(),
56 fPDG(0),
57 fCFManager(0x0),
58 fHistEventsProcessed(0x0),
a5d92cb4 59 fCorrelation(0x0),
f5ec6c30 60 fCountMC(0),
61 fCountAcc(0),
62 fCountReco(0),
63 fCountRecoAcc(0),
64 fCountRecoITSClusters(0),
65 fCountRecoPPR(0),
66 fEvents(0),
67 fFillFromGenerated(kFALSE),
a5d92cb4 68 fMinITSClusters(5),
69 fAcceptanceUnf(kTRUE)
f5ec6c30 70{
71 //
72 //Default ctor
73 //
74}
75//___________________________________________________________________________
76AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
77 AliAnalysisTaskSE(name),
78 fPDG(0),
79 fCFManager(0x0),
80 fHistEventsProcessed(0x0),
a5d92cb4 81 fCorrelation(0x0),
f5ec6c30 82 fCountMC(0),
83 fCountAcc(0),
84 fCountReco(0),
85 fCountRecoAcc(0),
86 fCountRecoITSClusters(0),
87 fCountRecoPPR(0),
88 fEvents(0),
89 fFillFromGenerated(kFALSE),
a5d92cb4 90 fMinITSClusters(5),
91 fAcceptanceUnf(kTRUE)
f5ec6c30 92{
93 //
94 // Constructor. Initialization of Inputs and Outputs
95 //
96 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
97 /*
98 DefineInput(0) and DefineOutput(0)
99 are taken care of by AliAnalysisTaskSE constructor
100 */
101 DefineOutput(1,TH1I::Class());
102 DefineOutput(2,AliCFContainer::Class());
a5d92cb4 103 DefineOutput(3,THnSparseD::Class());
f5ec6c30 104}
105
106//___________________________________________________________________________
107AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
108{
109 //
110 // Assignment operator
111 //
112 if (this!=&c) {
113 AliAnalysisTaskSE::operator=(c) ;
114 fPDG = c.fPDG;
115 fCFManager = c.fCFManager;
116 fHistEventsProcessed = c.fHistEventsProcessed;
117 }
118 return *this;
119}
120
121//___________________________________________________________________________
122AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
123 AliAnalysisTaskSE(c),
124 fPDG(c.fPDG),
125 fCFManager(c.fCFManager),
126 fHistEventsProcessed(c.fHistEventsProcessed),
a5d92cb4 127 fCorrelation(c.fCorrelation),
f5ec6c30 128 fCountMC(c.fCountMC),
129 fCountAcc(c.fCountAcc),
130 fCountReco(c.fCountReco),
131 fCountRecoAcc(c.fCountRecoAcc),
132 fCountRecoITSClusters(c.fCountRecoITSClusters),
133 fCountRecoPPR(c.fCountRecoPPR),
134 fEvents(c.fEvents),
135 fFillFromGenerated(c.fFillFromGenerated),
a5d92cb4 136 fMinITSClusters(c.fMinITSClusters),
137 fAcceptanceUnf(c.fAcceptanceUnf)
f5ec6c30 138{
139 //
140 // Copy Constructor
141 //
142}
143
144//___________________________________________________________________________
145AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
146 //
147 //destructor
148 //
149 if (fCFManager) delete fCFManager ;
150 if (fHistEventsProcessed) delete fHistEventsProcessed ;
a5d92cb4 151 if (fCorrelation) delete fCorrelation ;
f5ec6c30 152}
153
154//_________________________________________________
155void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
156{
157 //
158 // Main loop function
159 //
160
161 if (fFillFromGenerated){
162 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
163 }
164
165 if (!fInputEvent) {
166 Error("UserExec","NO EVENT FOUND!");
167 return;
168 }
169
170 fEvents++;
171 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
172 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
173 fCFManager->SetEventInfo(aodEvent);
174
175 // MC-event selection
176 Double_t containerInput[11] ;
177
178 //loop on the MC event
179
180 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
181 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
182 Int_t icountMC = 0;
183 Int_t icountAcc = 0;
184 Int_t icountReco = 0;
185 Int_t icountRecoAcc = 0;
186 Int_t icountRecoITSClusters = 0;
187 Int_t icountRecoPPR = 0;
188
e8158e56 189 Int_t cquarks = 0;
190
f5ec6c30 191 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
192 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
e8158e56 193 if (mcPart->GetPdgCode() == 4) cquarks++;
194 if (mcPart->GetPdgCode() == -4) cquarks++;
f5ec6c30 195 if (!mcPart) {
196 AliWarning("Particle not found in tree, skipping");
197 continue;
198 }
199
200 // check the MC-level cuts
e8158e56 201
f5ec6c30 202 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
e8158e56 203 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
204 Int_t abspdgGranma = TMath::Abs(pdgGranma);
205 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
206 AliInfo(Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
207 continue; // skipping particles that don't come from c quark
208 }
f5ec6c30 209
e8158e56 210 // if (TMath::Abs(pdgGranma)!=4) {
211
f5ec6c30 212 // fill the container for Gen-level selection
e8158e56 213 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
f5ec6c30 214 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
215 containerInput[0] = vectorMC[0];
216 containerInput[1] = vectorMC[1] ;
217 containerInput[2] = vectorMC[2] ;
218 containerInput[3] = vectorMC[3] ;
219 containerInput[4] = vectorMC[4] ;
220 containerInput[5] = vectorMC[5] ; // in micron
221 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
222 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
223 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
224 containerInput[9] = -100000.; // dummy value, meaningless in MC, in micron^2
225 containerInput[10] = 1.01; // dummy value, meaningless in MC
e8158e56 226 containerInput[11] = vectorMC[6]; // dummy value, meaningless in MC
f5ec6c30 227 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
228 icountMC++;
229
230 // check the MC-Acceptance level cuts
231 // since standard CF functions are not applicable, using Kine Cuts on daughters
232
233 Int_t daughter0 = mcPart->GetDaughter(0);
234 Int_t daughter1 = mcPart->GetDaughter(1);
235 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
236 if (daughter0 == 0 || daughter1 == 0) {
237 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
238 }
421623bd 239 if (TMath::Abs(daughter1 - daughter0) != 1) {
f5ec6c30 240 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
241 }
242 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
243 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
244 if (!mcPartDaughter0 || !mcPartDaughter1) {
245 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
246 }
247 Double_t eta0 = mcPartDaughter0->Eta();
248 Double_t eta1 = mcPartDaughter1->Eta();
249 Double_t y0 = mcPartDaughter0->Y();
250 Double_t y1 = mcPartDaughter1->Y();
251 Double_t pt0 = mcPartDaughter0->Pt();
252 Double_t pt1 = mcPartDaughter1->Pt();
253 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
254 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
255 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
256 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
257 if (daught0inAcceptance && daught1inAcceptance) {
258 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
259 AliDebug(2, "Daughter particles in acceptance");
260 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
261 AliInfo("Inconsistency with CF cut for daughter 0!");
262 }
263 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
264 AliInfo("Inconsistency with CF cut for daughter 1!");
265 }
266 fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptance);
267 icountAcc++;
268 }
269 }
270 else {
271 AliDebug(3,"Problems in filling the container");
272 continue;
273 }
e8158e56 274 }
275 if (cquarks<2) AliInfo(Form("Event found with %d c-quarks", cquarks));
f5ec6c30 276
277 AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC));
278 AliDebug(2, Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
279
280 // Now go to rec level
281 fCountMC += icountMC;
282 fCountAcc += icountAcc;
283
1e090d47 284 // AOD primary vertex
285 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
286
f5ec6c30 287 // load heavy flavour vertices
1e090d47 288 TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
289 if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
290 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
f5ec6c30 291
1e090d47 292 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
f5ec6c30 293
1e090d47 294 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
295 Bool_t unsetvtx=kFALSE;
296 if(!d0tokpi->GetOwnPrimaryVtx()) {
297 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
298 unsetvtx=kTRUE;
299 }
300
f5ec6c30 301 // find associated MC particle
1e090d47 302 Int_t mcLabel = d0tokpi->MatchToMC(421, mcArray) ;
f5ec6c30 303 if (mcLabel == -1)
304 {
305 AliDebug(2,"No MC particle found");
306 }
307 else {
308 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
309 if (!mcVtxHF) {
310 AliWarning("Could not find associated MC in AOD MC tree");
311 continue;
312 }
313
aed46c56 314 // check whether the daughters have kTPCrefit set
315 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
316 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
317 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit))) {
318 // skipping if at least one daughter does not have kTPCrefit
319 continue;
320 }
321
8855c601 322 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
aed46c56 323 Int_t okD0, okD0bar;
8855c601 324 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
aed46c56 325 // skipping candidate
326 continue;
327 }
328
f5ec6c30 329 // check if associated MC v0 passes the cuts
330 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
331 AliDebug(2, "Skipping the particles due to cuts");
332 continue;
333 }
e8158e56 334 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
335 Int_t abspdgGranma = TMath::Abs(pdgGranma);
336 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
337 AliInfo(Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
338 continue; // skipping particles that don't come from c quark
339 }
f5ec6c30 340
341 // fill the container...
1e090d47 342 Double_t pt = d0tokpi->Pt();
343 Double_t rapidity = d0tokpi->YD0();
f5ec6c30 344
345 Double_t cosThetaStar = 9999.;
346 Double_t pTpi = 0.;
347 Double_t pTK = 0.;
1e090d47 348 Double_t dca = d0tokpi->GetDCA();
f5ec6c30 349 Double_t d0pi = 0.;
350 Double_t d0K = 0.;
1e090d47 351 Double_t d0xd0 = d0tokpi->Prodd0d0();
e8158e56 352 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
353 Double_t phi = d0tokpi->Phi();
f5ec6c30 354 Int_t pdgCode = mcVtxHF->GetPdgCode();
355 if (pdgCode > 0){
1e090d47 356 cosThetaStar = d0tokpi->CosThetaStarD0();
357 pTpi = d0tokpi->PtProng(0);
358 pTK = d0tokpi->PtProng(1);
359 d0pi = d0tokpi->Getd0Prong(0);
360 d0K = d0tokpi->Getd0Prong(1);
f5ec6c30 361 }
362 else {
1e090d47 363 cosThetaStar = d0tokpi->CosThetaStarD0bar();
364 pTpi = d0tokpi->PtProng(1);
365 pTK = d0tokpi->PtProng(0);
366 d0pi = d0tokpi->Getd0Prong(1);
367 d0K = d0tokpi->Getd0Prong(0);
f5ec6c30 368 }
369
1e090d47 370 Double_t cT = d0tokpi->CtD0();
f5ec6c30 371
372 if (!fFillFromGenerated){
373 // ...either with reconstructed values....
374 containerInput[0] = pt;
375 containerInput[1] = rapidity;
376 containerInput[2] = cosThetaStar;
377 containerInput[3] = pTpi;
378 containerInput[4] = pTK;
379 containerInput[5] = cT*1.E4; // in micron
380 containerInput[6] = dca*1.E4; // in micron
381 containerInput[7] = d0pi*1.E4; // in micron
382 containerInput[8] = d0K*1.E4; // in micron
383 containerInput[9] = d0xd0*1.E8; // in micron^2
384 containerInput[10] = cosPointingAngle; // in micron
e8158e56 385 containerInput[11] = phi;
f5ec6c30 386 }
387 else {
388 // ... or with generated values
e8158e56 389 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
f5ec6c30 390 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
391 containerInput[0] = vectorMC[0];
392 containerInput[1] = vectorMC[1] ;
393 containerInput[2] = vectorMC[2] ;
394 containerInput[3] = vectorMC[3] ;
395 containerInput[4] = vectorMC[4] ;
396 containerInput[5] = vectorMC[5] ; // in micron
397 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
398 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
399 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
400 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
401 containerInput[10] = 1.01; // dummy value, meaningless in MC
e8158e56 402 containerInput[11] = vectorMC[6];
f5ec6c30 403 }
404 else {
405 AliDebug(3,"Problems in filling the container");
406 continue;
407 }
408 }
409 AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
410 icountReco++;
411 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
412
413 // cut in acceptance
1e090d47 414 Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
415 Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
f5ec6c30 416 if (acceptanceProng0 && acceptanceProng1) {
417 AliDebug(2,"D0 reco daughters in acceptance");
418 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
a5d92cb4 419 icountRecoAcc++;
420
421 if(fAcceptanceUnf){
422
423 Double_t fill[4]; //fill response matrix
424
425 // dimensions 0&1 : pt,eta (Rec)
426
427 fill[0] = pt ;
428 fill[1] = rapidity;
429
430 // dimensions 2&3 : pt,eta (MC)
431
432 fill[2] = mcVtxHF->Pt();
433 fill[3] = mcVtxHF->Y();
434
435 fCorrelation->Fill(fill);
436
437 }
f5ec6c30 438
439 // cut on the min n. of clusters in ITS
440 Int_t ncls0=0;
1e090d47 441 for(Int_t l=0;l<6;l++) if(TESTBIT(d0tokpi->GetITSClusterMap(),l)) ncls0++;
f5ec6c30 442 AliDebug(2, Form("n clusters = %d", ncls0));
443 if (ncls0 >= fMinITSClusters){
444 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
445 icountRecoITSClusters++;
446 AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
447
448 // PPR cuts
844d235c 449 Double_t cuts[6] = {9999999., 1.1, 0., 9999999., 9999999., 0.};
f5ec6c30 450 if (pt <= 1){
451 cuts[0] = 400;
452 cuts[1] = 0.8;
453 cuts[2] = 0.5;
454 cuts[3] = 500;
455 cuts[4] = -20000;
456 cuts[5] = 0.5;
457 }
458 else if (pt > 1 && pt <= 2){
459 cuts[0] = 300;
460 cuts[1] = 0.8;
461 cuts[2] = 0.6;
462 cuts[3] = 500;
463 cuts[4] = -20000;
464 cuts[5] = 0.6;
465 }
466 else if (pt > 2 && pt <= 3){
467 cuts[0] = 200;
468 cuts[1] = 0.8;
469 cuts[2] = 0.7;
470 cuts[3] = 500;
471 cuts[4] = -20000;
472 cuts[5] = 0.8;
473 }
474 else if (pt > 3 && pt <= 5){
475 cuts[0] = 200;
476 cuts[1] = 0.8;
477 cuts[2] = 0.7;
478 cuts[3] = 500;
479 cuts[4] = -10000;
480 cuts[5] = 0.8;
481 }
482 else if (pt > 5){
483 cuts[0] = 200;
484 cuts[1] = 0.8;
485 cuts[2] = 0.7;
486 cuts[3] = 500;
487 cuts[4] = -5000;
488 cuts[5] = 0.8;
489 }
490 if (dca*1E4 < cuts[0]
491 && TMath::Abs(cosThetaStar) < cuts[1]
492 && pTpi > cuts[2]
493 && pTK > cuts[2]
494 && TMath::Abs(d0pi*1E4) < cuts[3]
495 && TMath::Abs(d0K*1E4) < cuts[3]
496 && d0xd0*1E8 < cuts[4]
497 && cosPointingAngle > cuts[5]
498 ){
499
500 AliDebug(2,"Particle passed PPR cuts");
501 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
502 icountRecoPPR++;
a5d92cb4 503
504 if(!fAcceptanceUnf){ // unfolding
505
506 Double_t fill[4]; //fill response matrix
507
508 // dimensions 0&1 : pt,eta (Rec)
509
510 fill[0] = pt ;
511 fill[1] = rapidity;
512
513 // dimensions 2&3 : pt,eta (MC)
514
515 fill[2] = mcVtxHF->Pt();
516 fill[3] = mcVtxHF->Y();
517
518 fCorrelation->Fill(fill);
519
520 }
f5ec6c30 521 }
522 else{
523 AliDebug(2,"Particle skipped due to PPR cuts");
524 if (dca*1E4 > cuts[0]){
525 AliDebug(2,"Problems with dca");
526 }
527 if ( TMath::Abs(cosThetaStar) > cuts[1]){
528 AliDebug(2,"Problems with cosThetaStar");
529 }
530 if (pTpi < cuts[2]){
531 AliDebug(2,"Problems with pTpi");
532 }
533 if (pTK < cuts[2]){
534 AliDebug(2,"Problems with pTK");
535 }
536 if (TMath::Abs(d0pi*1E4) > cuts[3]){
537 AliDebug(2,"Problems with d0pi");
538 }
539 if (TMath::Abs(d0K*1E4) > cuts[3]){
540 AliDebug(2,"Problems with d0K");
541 }
542 if (d0xd0*1E8 > cuts[4]){
543 AliDebug(2,"Problems with d0xd0");
544 }
545 if (cosPointingAngle < cuts[5]){
546 AliDebug(2,"Problems with cosPointingAngle");
547 }
548 }
549 }
550 }
551 }
1e090d47 552 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
553 } // end loop on D0->Kpi
f5ec6c30 554
555 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
556
557 fCountReco+= icountReco;
558 fCountRecoAcc+= icountRecoAcc;
559 fCountRecoITSClusters+= icountRecoITSClusters;
560 fCountRecoPPR+= icountRecoPPR;
561
562 fHistEventsProcessed->Fill(0);
563 /* PostData(0) is taken care of by AliAnalysisTaskSE */
564 PostData(1,fHistEventsProcessed) ;
565 PostData(2,fCFManager->GetParticleContainer()) ;
a5d92cb4 566 PostData(3,fCorrelation) ;
f5ec6c30 567}
568
569
570//___________________________________________________________________________
571void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
572{
573 // The Terminate() function is the last function to be called during
574 // a query. It always runs on the client, it can be used to present
575 // the results graphically or save the results to file.
576
f5ec6c30 577 AliAnalysisTaskSE::Terminate();
578
579 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
580 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
581 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
582 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
583 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
584 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
585
586 // draw some example plots....
587
588 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
33aa234e 589 if(!cont) {
590 printf("CONTAINER NOT FOUND\n");
591 return;
592 }
f5ec6c30 593 // projecting the containers to obtain histograms
594 // first argument = variable, second argument = step
595
596 // MC-level
597 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
598 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
599 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
600 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
601 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
602 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
603 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
604 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
605 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
606 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
607 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
e8158e56 608 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
f5ec6c30 609
610 // MC-Acceptance level
611 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
612 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
613 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
614 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
615 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
616 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
617 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
618 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
619 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
620 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
621 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
e8158e56 622 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
f5ec6c30 623
624 // Reco-level
625 TH1D* h02 = cont->ShowProjection(0,2) ; // pt
626 TH1D* h12 = cont->ShowProjection(1,2) ; // rapidity
627 TH1D* h22 = cont->ShowProjection(2,2) ; // cosThetaStar
628 TH1D* h32 = cont->ShowProjection(3,2) ; // pTpi
629 TH1D* h42 = cont->ShowProjection(4,2) ; // pTK
630 TH1D* h52 = cont->ShowProjection(5,2) ; // cT
631 TH1D* h62 = cont->ShowProjection(6,2) ; // dca
632 TH1D* h72 = cont->ShowProjection(7,2) ; // d0pi
633 TH1D* h82 = cont->ShowProjection(8,2) ; // d0K
634 TH1D* h92 = cont->ShowProjection(9,2) ; // d0xd0
635 TH1D* h102 = cont->ShowProjection(10,2) ; // cosPointingAngle
e8158e56 636 TH1D* h112 = cont->ShowProjection(11,2) ; // phi
f5ec6c30 637
638 h00->SetTitle("pT_D0 (GeV/c)");
639 h10->SetTitle("rapidity");
640 h20->SetTitle("cosThetaStar");
641 h30->SetTitle("pT_pi (GeV/c)");
642 h40->SetTitle("pT_K (Gev/c)");
643 h50->SetTitle("cT (#mum)");
644 h60->SetTitle("dca (#mum)");
645 h70->SetTitle("d0_pi (#mum)");
646 h80->SetTitle("d0_K (#mum)");
647 h90->SetTitle("d0xd0 (#mum^2)");
648 h100->SetTitle("cosPointingAngle");
e8158e56 649 h100->SetTitle("phi (rad)");
f5ec6c30 650
651 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
652 h10->GetXaxis()->SetTitle("rapidity");
653 h20->GetXaxis()->SetTitle("cosThetaStar");
654 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
655 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
656 h50->GetXaxis()->SetTitle("cT (#mum)");
657 h60->GetXaxis()->SetTitle("dca (#mum)");
658 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
659 h80->GetXaxis()->SetTitle("d0_K (#mum)");
660 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
661 h100->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 662 h110->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 663
664 h01->SetTitle("pT_D0 (GeV/c)");
665 h11->SetTitle("rapidity");
666 h21->SetTitle("cosThetaStar");
667 h31->SetTitle("pT_pi (GeV/c)");
668 h41->SetTitle("pT_K (Gev/c)");
669 h51->SetTitle("cT (#mum)");
670 h61->SetTitle("dca (#mum)");
671 h71->SetTitle("d0_pi (#mum)");
672 h81->SetTitle("d0_K (#mum)");
673 h91->SetTitle("d0xd0 (#mum^2)");
674 h101->SetTitle("cosPointingAngle");
e8158e56 675 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 676
677 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
678 h11->GetXaxis()->SetTitle("rapidity");
679 h21->GetXaxis()->SetTitle("cosThetaStar");
680 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
681 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
682 h51->GetXaxis()->SetTitle("cT (#mum)");
683 h61->GetXaxis()->SetTitle("dca (#mum)");
684 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
685 h81->GetXaxis()->SetTitle("d0_K (#mum)");
686 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
687 h101->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 688 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 689
690 h02->SetTitle("pT_D0 (GeV/c)");
691 h12->SetTitle("rapidity");
692 h22->SetTitle("cosThetaStar");
693 h32->SetTitle("pT_pi (GeV/c)");
694 h42->SetTitle("pT_K (Gev/c)");
695 h52->SetTitle("cT (#mum)");
696 h62->SetTitle("dca (#mum)");
697 h72->SetTitle("d0_pi (#mum)");
698 h82->SetTitle("d0_K (#mum)");
699 h92->SetTitle("d0xd0 (#mum^2)");
700 h102->SetTitle("cosPointingAngle");
e8158e56 701 h112->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 702
703 h02->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
704 h12->GetXaxis()->SetTitle("rapidity");
705 h22->GetXaxis()->SetTitle("cosThetaStar");
706 h32->GetXaxis()->SetTitle("pT_pi (GeV/c)");
707 h42->GetXaxis()->SetTitle("pT_K (Gev/c)");
708 h52->GetXaxis()->SetTitle("cT (#mum)");
709 h62->GetXaxis()->SetTitle("dca (#mum)");
710 h72->GetXaxis()->SetTitle("d0_pi (#mum)");
711 h82->GetXaxis()->SetTitle("d0_K (#mum)");
712 h92->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
713 h102->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 714 h112->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 715
716 Double_t max0 = h00->GetMaximum();
717 Double_t max1 = h10->GetMaximum();
718 Double_t max2 = h20->GetMaximum();
719 Double_t max3 = h30->GetMaximum();
720 Double_t max4 = h40->GetMaximum();
721 Double_t max5 = h50->GetMaximum();
722 Double_t max6 = h60->GetMaximum();
723 Double_t max7 = h70->GetMaximum();
724 Double_t max8 = h80->GetMaximum();
725 Double_t max9 = h90->GetMaximum();
726 Double_t max10 = h100->GetMaximum();
e8158e56 727 Double_t max11 = h110->GetMaximum();
f5ec6c30 728
729 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
730 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
731 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
732 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
733 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
734 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
735 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
736 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
737 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
738 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
739 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 740 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 741
742 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
743 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
744 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
745 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
746 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
747 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
748 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
749 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
750 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
751 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
752 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 753 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 754
755 h00->SetMarkerStyle(20);
756 h10->SetMarkerStyle(24);
757 h20->SetMarkerStyle(21);
758 h30->SetMarkerStyle(25);
759 h40->SetMarkerStyle(27);
760 h50->SetMarkerStyle(28);
761 h60->SetMarkerStyle(20);
762 h70->SetMarkerStyle(24);
763 h80->SetMarkerStyle(21);
764 h90->SetMarkerStyle(25);
765 h100->SetMarkerStyle(27);
e8158e56 766 h110->SetMarkerStyle(28);
f5ec6c30 767
768 h00->SetMarkerColor(2);
769 h10->SetMarkerColor(2);
770 h20->SetMarkerColor(2);
771 h30->SetMarkerColor(2);
772 h40->SetMarkerColor(2);
773 h50->SetMarkerColor(2);
774 h60->SetMarkerColor(2);
775 h70->SetMarkerColor(2);
776 h80->SetMarkerColor(2);
777 h90->SetMarkerColor(2);
778 h100->SetMarkerColor(2);
e8158e56 779 h110->SetMarkerColor(2);
f5ec6c30 780
781 h01->SetMarkerStyle(20) ;
782 h11->SetMarkerStyle(24) ;
783 h21->SetMarkerStyle(21) ;
784 h31->SetMarkerStyle(25) ;
785 h41->SetMarkerStyle(27) ;
786 h51->SetMarkerStyle(28) ;
787 h61->SetMarkerStyle(20);
788 h71->SetMarkerStyle(24);
789 h81->SetMarkerStyle(21);
790 h91->SetMarkerStyle(25);
791 h101->SetMarkerStyle(27);
e8158e56 792 h111->SetMarkerStyle(28);
f5ec6c30 793
794 h01->SetMarkerColor(8);
795 h11->SetMarkerColor(8);
796 h21->SetMarkerColor(8);
797 h31->SetMarkerColor(8);
798 h41->SetMarkerColor(8);
799 h51->SetMarkerColor(8);
800 h61->SetMarkerColor(8);
801 h71->SetMarkerColor(8);
802 h81->SetMarkerColor(8);
803 h91->SetMarkerColor(8);
804 h101->SetMarkerColor(8);
e8158e56 805 h111->SetMarkerColor(8);
f5ec6c30 806
807 h02->SetMarkerStyle(20) ;
808 h12->SetMarkerStyle(24) ;
809 h22->SetMarkerStyle(21) ;
810 h32->SetMarkerStyle(25) ;
811 h42->SetMarkerStyle(27) ;
812 h52->SetMarkerStyle(28) ;
813 h62->SetMarkerStyle(20);
814 h72->SetMarkerStyle(24);
815 h82->SetMarkerStyle(21);
816 h92->SetMarkerStyle(25);
817 h102->SetMarkerStyle(27);
e8158e56 818 h112->SetMarkerStyle(28);
f5ec6c30 819
820 h02->SetMarkerColor(4);
821 h12->SetMarkerColor(4);
822 h22->SetMarkerColor(4);
823 h32->SetMarkerColor(4);
824 h42->SetMarkerColor(4);
825 h52->SetMarkerColor(4);
826 h62->SetMarkerColor(4);
827 h72->SetMarkerColor(4);
828 h82->SetMarkerColor(4);
829 h92->SetMarkerColor(4);
830 h102->SetMarkerColor(4);
e8158e56 831 h112->SetMarkerColor(4);
f5ec6c30 832
833 gStyle->SetCanvasColor(0);
834 gStyle->SetFrameFillColor(0);
835 gStyle->SetTitleFillColor(0);
836 gStyle->SetStatColor(0);
837
838 // drawing in 2 separate canvas for a matter of clearity
839 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
840 c1->Divide(3,3);
841
842 c1->cd(1);
843 h00->Draw("p");
844 c1->cd(1);
845 c1->cd(2);
846 h01->Draw("p");
847 c1->cd(2);
848 c1->cd(3);
849 h02->Draw("p");
850 c1->cd(3);
851 c1->cd(4);
852 h10->Draw("p");
853 c1->cd(4);
854 c1->cd(5);
855 h11->Draw("p");
856 c1->cd(5);
857 c1->cd(6);
858 h12->Draw("p");
859 c1->cd(6);
860 c1->cd(7);
861 h20->Draw("p");
862 c1->cd(7);
863 c1->cd(8);
864 h21->Draw("p");
865 c1->cd(8);
866 c1->cd(9);
867 h22->Draw("p");
868 c1->cd(9);
869 c1->cd();
870
871 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
872 c2->Divide(3,3);
873 c2->cd(1);
874 h30->Draw("p");
875 c2->cd(1);
876 c2->cd(2);
877 h31->Draw("p");
878 c2->cd(2);
879 c2->cd(3);
880 h32->Draw("p");
881 c2->cd(3);
882 c2->cd(4);
883 h40->Draw("p");
884 c2->cd(4);
885 c2->cd(5);
886 h41->Draw("p");
887 c2->cd(5);
888 c2->cd(6);
889 h42->Draw("p");
890 c2->cd(6);
891 c2->cd(7);
892 h50->Draw("p");
893 c2->cd(7);
894 c2->cd(8);
895 h51->Draw("p");
896 c2->cd(8);
897 c2->cd(9);
898 h52->Draw("p");
899 c2->cd(9);
900 c2->cd();
901
902 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
903 c3->Divide(3,3);
904 c3->cd(1);
905 h60->Draw("p");
906 c3->cd(1);
907 c3->cd(2);
908 h61->Draw("p");
909 c3->cd(2);
910 c3->cd(3);
911 h62->Draw("p");
912 c3->cd(3);
913 c3->cd(4);
914 h70->Draw("p");
915 c3->cd(4);
916 c3->cd(5);
917 h71->Draw("p");
918 c3->cd(5);
919 c3->cd(6);
920 h72->Draw("p");
921 c3->cd(6);
922 c3->cd(7);
923 h80->Draw("p");
924 c3->cd(7);
925 c3->cd(8);
926 h81->Draw("p");
927 c3->cd(8);
928 c3->cd(9);
929 h82->Draw("p");
930 c3->cd(9);
931 c3->cd();
932
e8158e56 933 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
934 c4->Divide(3,3);
f5ec6c30 935 c4->cd(1);
936 h90->Draw("p");
937 c4->cd(1);
938 c4->cd(2);
939 h91->Draw("p");
940 c4->cd(2);
941 c4->cd(3);
942 h92->Draw("p");
943 c4->cd(3);
944 c4->cd(4);
945 h100->Draw("p");
946 c4->cd(4);
947 c4->cd(5);
948 h101->Draw("p");
949 c4->cd(5);
950 c4->cd(6);
951 h102->Draw("p");
952 c4->cd(6);
e8158e56 953 c4->cd(7);
954 h110->Draw("p");
955 c4->cd(7);
956 c4->cd(8);
957 h111->Draw("p");
958 c4->cd(8);
959 c4->cd(9);
960 h112->Draw("p");
961 c4->cd(9);
f5ec6c30 962 c4->cd();
963
a5d92cb4 964 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
965
966 TH2D* corr1 =hcorr->Projection(0,2);
967 TH2D* corr2 = hcorr->Projection(1,3);
968
969 TCanvas * c7 =new TCanvas("c7","",800,400);
970 c7->Divide(2,1);
971 c7->cd(1);
972 corr1->Draw("text");
973 c7->cd(2);
974 corr2->Draw("text");
975
976
844d235c 977 TFile* file_projection = new TFile("file_projection.root","RECREATE");
a5d92cb4 978
979 corr1->Write();
980 corr2->Write();
844d235c 981 h00->Write("pT_D0_step0");
982 h10->Write("rapidity_step0");
983 h20->Write("cosThetaStar_step0");
984 h30->Write("pT_pi_step0");
985 h40->Write("pT_K_step0");
986 h50->Write("cT_step0");
987 h60->Write("dca_step0");
988 h70->Write("d0_pi_step0");
989 h80->Write("d0_K_step0");
990 h90->Write("d0xd0_step0");
991 h100->Write("cosPointingAngle_step0");
e8158e56 992 h110->Write("phi_step0");
844d235c 993
994 h01->Write("pT_D0_step1");
995 h11->Write("rapidity_step1");
996 h21->Write("cosThetaStar_step1");
997 h31->Write("pT_pi_step1");
998 h41->Write("pT_K_step1");
999 h51->Write("cT_step1");
1000 h61->Write("dca_step1");
1001 h71->Write("d0_pi_step1");
1002 h81->Write("d0_K_step1");
1003 h91->Write("d0xd0_step1");
1004 h101->Write("cosPointingAngle_step1");
e8158e56 1005 h111->Write("phi_step1");
844d235c 1006
1007 h02->Write("pT_D0_step2");
1008 h12->Write("rapidity_step2");
1009 h22->Write("cosThetaStar_step2");
1010 h32->Write("pT_pi_step2");
1011 h42->Write("pT_K_step2");
1012 h52->Write("cT_step2");
1013 h62->Write("dca_step2");
1014 h72->Write("d0_pi_step2");
1015 h80->Write("d0_K_step2");
1016 h92->Write("d0xd0_step2");
1017 h102->Write("cosPointingAngle_step2");
e8158e56 1018 h112->Write("phi_step2");
844d235c 1019
1020 file_projection->Close();
1021
a5d92cb4 1022
1023
f5ec6c30 1024 /*
1025 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1026 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1027 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1028 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1029
1030 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1031 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1032 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1033 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1034 */
1035}
1036
1037//___________________________________________________________________________
1038void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1039 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1040 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1041 //
1042 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1043
1044 //slot #1
1045 OpenFile(1);
1046 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
1047}
1048
1049//___________________________________________________________________________
1050Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1051
1052 //
1053 // to calculate cos(ThetaStar) for generated particle
1054 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1055 // (see where the function is called)
1056 //
1057
1058 Int_t pdgvtx = mcPart->GetPdgCode();
1059 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1060 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1061 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1062 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1063 AliDebug(2,"This is a D0");
1064 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1065 mcPartDaughter0 = mcPartDaughter1;
1066 mcPartDaughter1 = mcPartdummy;
1067 }
1068 else{
1069 AliInfo("D0bar");
1070 }
1071 */
1072 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1073 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1074 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1075 AliDebug(2,"D0");
1076 }
1077 else{
1078 AliDebug(2,"D0bar");
1079 }
1080 if (pdgprong0 == 211){
1081 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1082 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1083 mcPartDaughter0 = mcPartDaughter1;
1084 mcPartDaughter1 = mcPartdummy;
1085 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1086 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1087 }
1088
1089 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1090 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1091 Double_t massp[2];
1092 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1093 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1094
1095 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
1096
1097 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1098 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1099 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1100 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1101 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1102 Double_t beta = p/e;
1103 Double_t gamma = e/massvtx;
1104 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1105 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1106
1107 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1108
1109 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1110 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1111 AliDebug(2,Form("cts = %f", cts));
1112 return cts;
1113}
1114//___________________________________________________________________________
1115Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1116
1117 //
1118 // to calculate cT for generated particle
1119 //
1120
1121 Double_t xmcPart[3] = {0,0,0};
1122 Double_t xdaughter0[3] = {0,0,0};
1123 Double_t xdaughter1[3] = {0,0,0};
1124 mcPart->XvYvZv(xmcPart); // cm
1125 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1126 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1127 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1128 xmcPart[1]*xmcPart[1]+
1129 xmcPart[2]*xmcPart[2]);
1130 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1131 xdaughter0[1]*xdaughter0[1]+
1132 xdaughter0[2]*xdaughter0[2]);
1133 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1134 xdaughter1[1]*xdaughter1[1]+
1135 xdaughter1[2]*xdaughter1[2]);
1136
1137 AliDebug(2, Form("D0: x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4));
1138 AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4));
1139 AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4));
1140
1141 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1142 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1143 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1144
1145 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1146 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1147 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1148
e8158e56 1149 if ((cT0 - cT1)>1E-5) {
1150 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
f5ec6c30 1151 }
1152
1153 // calculating cT from cT0
1154
1155 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1156 Double_t p = mcPart-> P();
1157 Double_t cT = cT0*mass/p;
1158 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1159 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1160 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1161 return cT;
1162}
1163//________________________________________________________________________________
1164Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1165
1166 //
1167 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1168 //
1169
1170 Bool_t isOk = kFALSE;
1171
1172 // check whether the D0 decays in pi+K
1173 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1174 // could use a cut in the CF, but not clear how to define a TDecayChannel
1175 // implemented for the time being as a cut on the number of daughters - see later when
1176 // getting the daughters
1177
1178 // getting the daughters
1179 Int_t daughter0 = mcPart->GetDaughter(0);
1180 Int_t daughter1 = mcPart->GetDaughter(1);
1181 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1182 if (daughter0 == 0 || daughter1 == 0) {
1183 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1184 return isOk;
1185 }
421623bd 1186 if (TMath::Abs(daughter1 - daughter0) != 1) {
f5ec6c30 1187 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1188 return isOk;
1189 }
1190 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1191 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1192 if (!mcPartDaughter0 || !mcPartDaughter1) {
1193 AliWarning("At least one Daughter Particle not found in tree, skipping");
1194 return isOk;
1195 }
2bcc508a 1196 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1197 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1198 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1199 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1200 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1201 return isOk;
1202 }
1203
f5ec6c30 1204 Double_t vtx1[3] = {0,0,0}; // primary vertex
1205 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1206 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1207 mcPart->XvYvZv(vtx1); // cm
1208 // getting vertex from daughters
1209 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1210 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1211 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1212 AliError("Daughters have different secondary vertex, skipping the track");
1213 return isOk;
1214 }
1215 Int_t nprongs = 2;
1216 Short_t charge = 0;
1217 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1218 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1219 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1220 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1221 // inverting in case the positive daughter is the second one
1222 positiveDaugh = mcPartDaughter1;
1223 negativeDaugh = mcPartDaughter0;
1224 }
1225 // getting the momentum from the daughters
1226 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1227 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1228 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1229
1230 Double_t d0[2] = {0.,0.};
1231
1232 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1233
1234 Double_t cosThetaStar = 0.;
1235 Double_t cosThetaStarD0 = 0.;
1236 Double_t cosThetaStarD0bar = 0.;
1237 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1238 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1239 if (mcPart->GetPdgCode() == 421){ // D0
1240 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1241 cosThetaStar = cosThetaStarD0;
1242 }
1243 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1244 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1245 cosThetaStar = cosThetaStarD0bar;
1246 }
1247 else{
1248 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1249 return vectorMC;
1250 }
1251 if (cosThetaStar < -1 || cosThetaStar > 1) {
1252 AliWarning("Invalid value for cosine Theta star, returning");
1253 return isOk;
1254 }
1255
1256 // calculate cos(Theta*) according to the method implemented herein
1257
1258 Double_t cts = 9999.;
1259 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1260 if (cts < -1 || cts > 1) {
1261 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1262 }
1263 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1264 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1265 }
1266
1267 Double_t cT = decay->Ct(421);
1268
1269 // calculate cT from the method implemented herein
1270 Double_t cT1 = 0.;
1271 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1272
1273 if (TMath::Abs(cT1 - cT)>0.001) {
1274 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1275 }
1276
1277 // get the pT of the daughters
1278
1279 Double_t pTpi = 0.;
1280 Double_t pTK = 0.;
1281
1282 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1283 pTpi = mcPartDaughter0->Pt();
1284 pTK = mcPartDaughter1->Pt();
1285 }
1286 else {
1287 pTpi = mcPartDaughter1->Pt();
1288 pTK = mcPartDaughter0->Pt();
1289 }
1290
1291 vectorMC[0] = mcPart->Pt();
1292 vectorMC[1] = mcPart->Y() ;
1293 vectorMC[2] = cosThetaStar ;
1294 vectorMC[3] = pTpi ;
1295 vectorMC[4] = pTK ;
1296 vectorMC[5] = cT*1.E4 ; // in micron
e8158e56 1297 vectorMC[6] = mcPart->Phi() ;
f5ec6c30 1298 isOk = kTRUE;
1299 return isOk;
1300}
e8158e56 1301//_________________________________________________________________________________________________
1302Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1303
1304 //
1305 // checking whether the very mother of the D0 is a charm or a bottom quark
1306 //
1307
1308 Int_t pdgGranma = 0;
1309 Int_t mother = 0;
1310 mother = mcPart->GetMother();
1311 Int_t istep = 0;
1312 while (mother >0 ){
1313 istep++;
1314 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1315 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1316 pdgGranma = mcGranma->GetPdgCode();
1317 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1318 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1319 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1320 break;
1321 }
1322 mother = mcGranma->GetMother();
1323 }
1324 return pdgGranma;
1325}
f5ec6c30 1326