]> 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
Adding root flags to C flags
[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
feedec5e 20// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
f5ec6c30 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
feedec5e 176 Double_t containerInput[12] ;
f5ec6c30 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)) {
caf4ca8b 206 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
e8158e56 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 }
caf4ca8b 275 if (cquarks<2) AliDebug(2,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 }
caf4ca8b 313 // check whether the daughters are K- and pi+
314 AliAODMCParticle* dg0MC=(AliAODMCParticle*)mcArray->At(mcVtxHF->GetDaughter(0));
315 AliAODMCParticle* dg1MC=(AliAODMCParticle*)mcArray->At(mcVtxHF->GetDaughter(1));
316 if(!(TMath::Abs(dg0MC->GetPdgCode())==321 && TMath::Abs(dg1MC->GetPdgCode())==211) &&
317 !(TMath::Abs(dg0MC->GetPdgCode())==211 && TMath::Abs(dg1MC->GetPdgCode())==321)) continue;
318
aed46c56 319 // check whether the daughters have kTPCrefit set
320 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
321 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
322 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit))) {
323 // skipping if at least one daughter does not have kTPCrefit
324 continue;
325 }
326
8855c601 327 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
aed46c56 328 Int_t okD0, okD0bar;
8855c601 329 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
aed46c56 330 // skipping candidate
331 continue;
332 }
333
f5ec6c30 334 // check if associated MC v0 passes the cuts
335 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
336 AliDebug(2, "Skipping the particles due to cuts");
337 continue;
338 }
e8158e56 339 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
340 Int_t abspdgGranma = TMath::Abs(pdgGranma);
341 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
caf4ca8b 342 AliDebug(2,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));
e8158e56 343 continue; // skipping particles that don't come from c quark
344 }
f5ec6c30 345
346 // fill the container...
1e090d47 347 Double_t pt = d0tokpi->Pt();
348 Double_t rapidity = d0tokpi->YD0();
f5ec6c30 349
350 Double_t cosThetaStar = 9999.;
351 Double_t pTpi = 0.;
352 Double_t pTK = 0.;
1e090d47 353 Double_t dca = d0tokpi->GetDCA();
f5ec6c30 354 Double_t d0pi = 0.;
355 Double_t d0K = 0.;
1e090d47 356 Double_t d0xd0 = d0tokpi->Prodd0d0();
e8158e56 357 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
358 Double_t phi = d0tokpi->Phi();
f5ec6c30 359 Int_t pdgCode = mcVtxHF->GetPdgCode();
360 if (pdgCode > 0){
1e090d47 361 cosThetaStar = d0tokpi->CosThetaStarD0();
362 pTpi = d0tokpi->PtProng(0);
363 pTK = d0tokpi->PtProng(1);
364 d0pi = d0tokpi->Getd0Prong(0);
365 d0K = d0tokpi->Getd0Prong(1);
f5ec6c30 366 }
367 else {
1e090d47 368 cosThetaStar = d0tokpi->CosThetaStarD0bar();
369 pTpi = d0tokpi->PtProng(1);
370 pTK = d0tokpi->PtProng(0);
371 d0pi = d0tokpi->Getd0Prong(1);
372 d0K = d0tokpi->Getd0Prong(0);
f5ec6c30 373 }
374
1e090d47 375 Double_t cT = d0tokpi->CtD0();
f5ec6c30 376
377 if (!fFillFromGenerated){
378 // ...either with reconstructed values....
379 containerInput[0] = pt;
380 containerInput[1] = rapidity;
381 containerInput[2] = cosThetaStar;
382 containerInput[3] = pTpi;
383 containerInput[4] = pTK;
384 containerInput[5] = cT*1.E4; // in micron
385 containerInput[6] = dca*1.E4; // in micron
386 containerInput[7] = d0pi*1.E4; // in micron
387 containerInput[8] = d0K*1.E4; // in micron
388 containerInput[9] = d0xd0*1.E8; // in micron^2
389 containerInput[10] = cosPointingAngle; // in micron
e8158e56 390 containerInput[11] = phi;
f5ec6c30 391 }
392 else {
393 // ... or with generated values
e8158e56 394 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
f5ec6c30 395 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
396 containerInput[0] = vectorMC[0];
397 containerInput[1] = vectorMC[1] ;
398 containerInput[2] = vectorMC[2] ;
399 containerInput[3] = vectorMC[3] ;
400 containerInput[4] = vectorMC[4] ;
401 containerInput[5] = vectorMC[5] ; // in micron
402 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
403 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
404 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
405 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
406 containerInput[10] = 1.01; // dummy value, meaningless in MC
e8158e56 407 containerInput[11] = vectorMC[6];
f5ec6c30 408 }
409 else {
410 AliDebug(3,"Problems in filling the container");
411 continue;
412 }
413 }
414 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]));
415 icountReco++;
416 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
417
418 // cut in acceptance
1e090d47 419 Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
420 Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
f5ec6c30 421 if (acceptanceProng0 && acceptanceProng1) {
422 AliDebug(2,"D0 reco daughters in acceptance");
423 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
a5d92cb4 424 icountRecoAcc++;
425
426 if(fAcceptanceUnf){
427
428 Double_t fill[4]; //fill response matrix
429
430 // dimensions 0&1 : pt,eta (Rec)
431
432 fill[0] = pt ;
433 fill[1] = rapidity;
434
435 // dimensions 2&3 : pt,eta (MC)
436
437 fill[2] = mcVtxHF->Pt();
438 fill[3] = mcVtxHF->Y();
439
440 fCorrelation->Fill(fill);
441
442 }
f5ec6c30 443
444 // cut on the min n. of clusters in ITS
445 Int_t ncls0=0;
1e090d47 446 for(Int_t l=0;l<6;l++) if(TESTBIT(d0tokpi->GetITSClusterMap(),l)) ncls0++;
f5ec6c30 447 AliDebug(2, Form("n clusters = %d", ncls0));
448 if (ncls0 >= fMinITSClusters){
449 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
450 icountRecoITSClusters++;
451 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));
452
453 // PPR cuts
844d235c 454 Double_t cuts[6] = {9999999., 1.1, 0., 9999999., 9999999., 0.};
f5ec6c30 455 if (pt <= 1){
456 cuts[0] = 400;
457 cuts[1] = 0.8;
458 cuts[2] = 0.5;
459 cuts[3] = 500;
460 cuts[4] = -20000;
461 cuts[5] = 0.5;
462 }
463 else if (pt > 1 && pt <= 2){
464 cuts[0] = 300;
465 cuts[1] = 0.8;
466 cuts[2] = 0.6;
467 cuts[3] = 500;
468 cuts[4] = -20000;
469 cuts[5] = 0.6;
470 }
471 else if (pt > 2 && pt <= 3){
472 cuts[0] = 200;
473 cuts[1] = 0.8;
474 cuts[2] = 0.7;
475 cuts[3] = 500;
476 cuts[4] = -20000;
477 cuts[5] = 0.8;
478 }
479 else if (pt > 3 && pt <= 5){
480 cuts[0] = 200;
481 cuts[1] = 0.8;
482 cuts[2] = 0.7;
483 cuts[3] = 500;
484 cuts[4] = -10000;
485 cuts[5] = 0.8;
486 }
487 else if (pt > 5){
488 cuts[0] = 200;
489 cuts[1] = 0.8;
490 cuts[2] = 0.7;
491 cuts[3] = 500;
492 cuts[4] = -5000;
493 cuts[5] = 0.8;
494 }
495 if (dca*1E4 < cuts[0]
496 && TMath::Abs(cosThetaStar) < cuts[1]
497 && pTpi > cuts[2]
498 && pTK > cuts[2]
499 && TMath::Abs(d0pi*1E4) < cuts[3]
500 && TMath::Abs(d0K*1E4) < cuts[3]
501 && d0xd0*1E8 < cuts[4]
502 && cosPointingAngle > cuts[5]
503 ){
504
505 AliDebug(2,"Particle passed PPR cuts");
506 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
507 icountRecoPPR++;
a5d92cb4 508
509 if(!fAcceptanceUnf){ // unfolding
510
511 Double_t fill[4]; //fill response matrix
512
513 // dimensions 0&1 : pt,eta (Rec)
514
515 fill[0] = pt ;
516 fill[1] = rapidity;
517
518 // dimensions 2&3 : pt,eta (MC)
519
520 fill[2] = mcVtxHF->Pt();
521 fill[3] = mcVtxHF->Y();
522
523 fCorrelation->Fill(fill);
524
525 }
f5ec6c30 526 }
527 else{
528 AliDebug(2,"Particle skipped due to PPR cuts");
529 if (dca*1E4 > cuts[0]){
530 AliDebug(2,"Problems with dca");
531 }
532 if ( TMath::Abs(cosThetaStar) > cuts[1]){
533 AliDebug(2,"Problems with cosThetaStar");
534 }
535 if (pTpi < cuts[2]){
536 AliDebug(2,"Problems with pTpi");
537 }
538 if (pTK < cuts[2]){
539 AliDebug(2,"Problems with pTK");
540 }
541 if (TMath::Abs(d0pi*1E4) > cuts[3]){
542 AliDebug(2,"Problems with d0pi");
543 }
544 if (TMath::Abs(d0K*1E4) > cuts[3]){
545 AliDebug(2,"Problems with d0K");
546 }
547 if (d0xd0*1E8 > cuts[4]){
548 AliDebug(2,"Problems with d0xd0");
549 }
550 if (cosPointingAngle < cuts[5]){
551 AliDebug(2,"Problems with cosPointingAngle");
552 }
553 }
554 }
555 }
556 }
1e090d47 557 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
558 } // end loop on D0->Kpi
f5ec6c30 559
560 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
561
562 fCountReco+= icountReco;
563 fCountRecoAcc+= icountRecoAcc;
564 fCountRecoITSClusters+= icountRecoITSClusters;
565 fCountRecoPPR+= icountRecoPPR;
566
567 fHistEventsProcessed->Fill(0);
568 /* PostData(0) is taken care of by AliAnalysisTaskSE */
569 PostData(1,fHistEventsProcessed) ;
570 PostData(2,fCFManager->GetParticleContainer()) ;
a5d92cb4 571 PostData(3,fCorrelation) ;
f5ec6c30 572}
573
574
575//___________________________________________________________________________
576void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
577{
578 // The Terminate() function is the last function to be called during
579 // a query. It always runs on the client, it can be used to present
580 // the results graphically or save the results to file.
581
f5ec6c30 582 AliAnalysisTaskSE::Terminate();
583
584 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
585 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
586 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
587 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));
588 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));
589 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
590
591 // draw some example plots....
592
593 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
33aa234e 594 if(!cont) {
595 printf("CONTAINER NOT FOUND\n");
596 return;
597 }
f5ec6c30 598 // projecting the containers to obtain histograms
599 // first argument = variable, second argument = step
600
601 // MC-level
602 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
603 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
604 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
605 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
606 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
607 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
608 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
609 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
610 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
611 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
612 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
e8158e56 613 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
f5ec6c30 614
615 // MC-Acceptance level
616 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
617 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
618 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
619 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
620 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
621 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
622 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
623 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
624 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
625 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
626 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
e8158e56 627 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
f5ec6c30 628
629 // Reco-level
630 TH1D* h02 = cont->ShowProjection(0,2) ; // pt
631 TH1D* h12 = cont->ShowProjection(1,2) ; // rapidity
632 TH1D* h22 = cont->ShowProjection(2,2) ; // cosThetaStar
633 TH1D* h32 = cont->ShowProjection(3,2) ; // pTpi
634 TH1D* h42 = cont->ShowProjection(4,2) ; // pTK
635 TH1D* h52 = cont->ShowProjection(5,2) ; // cT
636 TH1D* h62 = cont->ShowProjection(6,2) ; // dca
637 TH1D* h72 = cont->ShowProjection(7,2) ; // d0pi
638 TH1D* h82 = cont->ShowProjection(8,2) ; // d0K
639 TH1D* h92 = cont->ShowProjection(9,2) ; // d0xd0
640 TH1D* h102 = cont->ShowProjection(10,2) ; // cosPointingAngle
e8158e56 641 TH1D* h112 = cont->ShowProjection(11,2) ; // phi
f5ec6c30 642
643 h00->SetTitle("pT_D0 (GeV/c)");
644 h10->SetTitle("rapidity");
645 h20->SetTitle("cosThetaStar");
646 h30->SetTitle("pT_pi (GeV/c)");
647 h40->SetTitle("pT_K (Gev/c)");
648 h50->SetTitle("cT (#mum)");
649 h60->SetTitle("dca (#mum)");
650 h70->SetTitle("d0_pi (#mum)");
651 h80->SetTitle("d0_K (#mum)");
652 h90->SetTitle("d0xd0 (#mum^2)");
653 h100->SetTitle("cosPointingAngle");
e8158e56 654 h100->SetTitle("phi (rad)");
f5ec6c30 655
656 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
657 h10->GetXaxis()->SetTitle("rapidity");
658 h20->GetXaxis()->SetTitle("cosThetaStar");
659 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
660 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
661 h50->GetXaxis()->SetTitle("cT (#mum)");
662 h60->GetXaxis()->SetTitle("dca (#mum)");
663 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
664 h80->GetXaxis()->SetTitle("d0_K (#mum)");
665 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
666 h100->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 667 h110->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 668
669 h01->SetTitle("pT_D0 (GeV/c)");
670 h11->SetTitle("rapidity");
671 h21->SetTitle("cosThetaStar");
672 h31->SetTitle("pT_pi (GeV/c)");
673 h41->SetTitle("pT_K (Gev/c)");
674 h51->SetTitle("cT (#mum)");
675 h61->SetTitle("dca (#mum)");
676 h71->SetTitle("d0_pi (#mum)");
677 h81->SetTitle("d0_K (#mum)");
678 h91->SetTitle("d0xd0 (#mum^2)");
679 h101->SetTitle("cosPointingAngle");
e8158e56 680 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 681
682 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
683 h11->GetXaxis()->SetTitle("rapidity");
684 h21->GetXaxis()->SetTitle("cosThetaStar");
685 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
686 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
687 h51->GetXaxis()->SetTitle("cT (#mum)");
688 h61->GetXaxis()->SetTitle("dca (#mum)");
689 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
690 h81->GetXaxis()->SetTitle("d0_K (#mum)");
691 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
692 h101->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 693 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 694
695 h02->SetTitle("pT_D0 (GeV/c)");
696 h12->SetTitle("rapidity");
697 h22->SetTitle("cosThetaStar");
698 h32->SetTitle("pT_pi (GeV/c)");
699 h42->SetTitle("pT_K (Gev/c)");
700 h52->SetTitle("cT (#mum)");
701 h62->SetTitle("dca (#mum)");
702 h72->SetTitle("d0_pi (#mum)");
703 h82->SetTitle("d0_K (#mum)");
704 h92->SetTitle("d0xd0 (#mum^2)");
705 h102->SetTitle("cosPointingAngle");
e8158e56 706 h112->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 707
708 h02->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
709 h12->GetXaxis()->SetTitle("rapidity");
710 h22->GetXaxis()->SetTitle("cosThetaStar");
711 h32->GetXaxis()->SetTitle("pT_pi (GeV/c)");
712 h42->GetXaxis()->SetTitle("pT_K (Gev/c)");
713 h52->GetXaxis()->SetTitle("cT (#mum)");
714 h62->GetXaxis()->SetTitle("dca (#mum)");
715 h72->GetXaxis()->SetTitle("d0_pi (#mum)");
716 h82->GetXaxis()->SetTitle("d0_K (#mum)");
717 h92->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
718 h102->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 719 h112->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 720
721 Double_t max0 = h00->GetMaximum();
722 Double_t max1 = h10->GetMaximum();
723 Double_t max2 = h20->GetMaximum();
724 Double_t max3 = h30->GetMaximum();
725 Double_t max4 = h40->GetMaximum();
726 Double_t max5 = h50->GetMaximum();
727 Double_t max6 = h60->GetMaximum();
728 Double_t max7 = h70->GetMaximum();
729 Double_t max8 = h80->GetMaximum();
730 Double_t max9 = h90->GetMaximum();
731 Double_t max10 = h100->GetMaximum();
e8158e56 732 Double_t max11 = h110->GetMaximum();
f5ec6c30 733
734 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
735 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
736 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
737 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
738 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
739 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
740 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
741 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
742 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
743 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
744 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 745 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 746
747 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
748 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
749 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
750 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
751 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
752 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
753 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
754 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
755 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
756 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
757 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 758 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 759
760 h00->SetMarkerStyle(20);
761 h10->SetMarkerStyle(24);
762 h20->SetMarkerStyle(21);
763 h30->SetMarkerStyle(25);
764 h40->SetMarkerStyle(27);
765 h50->SetMarkerStyle(28);
766 h60->SetMarkerStyle(20);
767 h70->SetMarkerStyle(24);
768 h80->SetMarkerStyle(21);
769 h90->SetMarkerStyle(25);
770 h100->SetMarkerStyle(27);
e8158e56 771 h110->SetMarkerStyle(28);
f5ec6c30 772
773 h00->SetMarkerColor(2);
774 h10->SetMarkerColor(2);
775 h20->SetMarkerColor(2);
776 h30->SetMarkerColor(2);
777 h40->SetMarkerColor(2);
778 h50->SetMarkerColor(2);
779 h60->SetMarkerColor(2);
780 h70->SetMarkerColor(2);
781 h80->SetMarkerColor(2);
782 h90->SetMarkerColor(2);
783 h100->SetMarkerColor(2);
e8158e56 784 h110->SetMarkerColor(2);
f5ec6c30 785
786 h01->SetMarkerStyle(20) ;
787 h11->SetMarkerStyle(24) ;
788 h21->SetMarkerStyle(21) ;
789 h31->SetMarkerStyle(25) ;
790 h41->SetMarkerStyle(27) ;
791 h51->SetMarkerStyle(28) ;
792 h61->SetMarkerStyle(20);
793 h71->SetMarkerStyle(24);
794 h81->SetMarkerStyle(21);
795 h91->SetMarkerStyle(25);
796 h101->SetMarkerStyle(27);
e8158e56 797 h111->SetMarkerStyle(28);
f5ec6c30 798
799 h01->SetMarkerColor(8);
800 h11->SetMarkerColor(8);
801 h21->SetMarkerColor(8);
802 h31->SetMarkerColor(8);
803 h41->SetMarkerColor(8);
804 h51->SetMarkerColor(8);
805 h61->SetMarkerColor(8);
806 h71->SetMarkerColor(8);
807 h81->SetMarkerColor(8);
808 h91->SetMarkerColor(8);
809 h101->SetMarkerColor(8);
e8158e56 810 h111->SetMarkerColor(8);
f5ec6c30 811
812 h02->SetMarkerStyle(20) ;
813 h12->SetMarkerStyle(24) ;
814 h22->SetMarkerStyle(21) ;
815 h32->SetMarkerStyle(25) ;
816 h42->SetMarkerStyle(27) ;
817 h52->SetMarkerStyle(28) ;
818 h62->SetMarkerStyle(20);
819 h72->SetMarkerStyle(24);
820 h82->SetMarkerStyle(21);
821 h92->SetMarkerStyle(25);
822 h102->SetMarkerStyle(27);
e8158e56 823 h112->SetMarkerStyle(28);
f5ec6c30 824
825 h02->SetMarkerColor(4);
826 h12->SetMarkerColor(4);
827 h22->SetMarkerColor(4);
828 h32->SetMarkerColor(4);
829 h42->SetMarkerColor(4);
830 h52->SetMarkerColor(4);
831 h62->SetMarkerColor(4);
832 h72->SetMarkerColor(4);
833 h82->SetMarkerColor(4);
834 h92->SetMarkerColor(4);
835 h102->SetMarkerColor(4);
e8158e56 836 h112->SetMarkerColor(4);
f5ec6c30 837
838 gStyle->SetCanvasColor(0);
839 gStyle->SetFrameFillColor(0);
840 gStyle->SetTitleFillColor(0);
841 gStyle->SetStatColor(0);
842
843 // drawing in 2 separate canvas for a matter of clearity
844 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
845 c1->Divide(3,3);
846
847 c1->cd(1);
848 h00->Draw("p");
849 c1->cd(1);
850 c1->cd(2);
851 h01->Draw("p");
852 c1->cd(2);
853 c1->cd(3);
854 h02->Draw("p");
855 c1->cd(3);
856 c1->cd(4);
857 h10->Draw("p");
858 c1->cd(4);
859 c1->cd(5);
860 h11->Draw("p");
861 c1->cd(5);
862 c1->cd(6);
863 h12->Draw("p");
864 c1->cd(6);
865 c1->cd(7);
866 h20->Draw("p");
867 c1->cd(7);
868 c1->cd(8);
869 h21->Draw("p");
870 c1->cd(8);
871 c1->cd(9);
872 h22->Draw("p");
873 c1->cd(9);
874 c1->cd();
875
876 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
877 c2->Divide(3,3);
878 c2->cd(1);
879 h30->Draw("p");
880 c2->cd(1);
881 c2->cd(2);
882 h31->Draw("p");
883 c2->cd(2);
884 c2->cd(3);
885 h32->Draw("p");
886 c2->cd(3);
887 c2->cd(4);
888 h40->Draw("p");
889 c2->cd(4);
890 c2->cd(5);
891 h41->Draw("p");
892 c2->cd(5);
893 c2->cd(6);
894 h42->Draw("p");
895 c2->cd(6);
896 c2->cd(7);
897 h50->Draw("p");
898 c2->cd(7);
899 c2->cd(8);
900 h51->Draw("p");
901 c2->cd(8);
902 c2->cd(9);
903 h52->Draw("p");
904 c2->cd(9);
905 c2->cd();
906
907 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
908 c3->Divide(3,3);
909 c3->cd(1);
910 h60->Draw("p");
911 c3->cd(1);
912 c3->cd(2);
913 h61->Draw("p");
914 c3->cd(2);
915 c3->cd(3);
916 h62->Draw("p");
917 c3->cd(3);
918 c3->cd(4);
919 h70->Draw("p");
920 c3->cd(4);
921 c3->cd(5);
922 h71->Draw("p");
923 c3->cd(5);
924 c3->cd(6);
925 h72->Draw("p");
926 c3->cd(6);
927 c3->cd(7);
928 h80->Draw("p");
929 c3->cd(7);
930 c3->cd(8);
931 h81->Draw("p");
932 c3->cd(8);
933 c3->cd(9);
934 h82->Draw("p");
935 c3->cd(9);
936 c3->cd();
937
e8158e56 938 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
939 c4->Divide(3,3);
f5ec6c30 940 c4->cd(1);
941 h90->Draw("p");
942 c4->cd(1);
943 c4->cd(2);
944 h91->Draw("p");
945 c4->cd(2);
946 c4->cd(3);
947 h92->Draw("p");
948 c4->cd(3);
949 c4->cd(4);
950 h100->Draw("p");
951 c4->cd(4);
952 c4->cd(5);
953 h101->Draw("p");
954 c4->cd(5);
955 c4->cd(6);
956 h102->Draw("p");
957 c4->cd(6);
e8158e56 958 c4->cd(7);
959 h110->Draw("p");
960 c4->cd(7);
961 c4->cd(8);
962 h111->Draw("p");
963 c4->cd(8);
964 c4->cd(9);
965 h112->Draw("p");
966 c4->cd(9);
f5ec6c30 967 c4->cd();
968
a5d92cb4 969 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
970
971 TH2D* corr1 =hcorr->Projection(0,2);
972 TH2D* corr2 = hcorr->Projection(1,3);
973
974 TCanvas * c7 =new TCanvas("c7","",800,400);
975 c7->Divide(2,1);
976 c7->cd(1);
977 corr1->Draw("text");
978 c7->cd(2);
979 corr2->Draw("text");
980
981
844d235c 982 TFile* file_projection = new TFile("file_projection.root","RECREATE");
a5d92cb4 983
984 corr1->Write();
985 corr2->Write();
844d235c 986 h00->Write("pT_D0_step0");
987 h10->Write("rapidity_step0");
988 h20->Write("cosThetaStar_step0");
989 h30->Write("pT_pi_step0");
990 h40->Write("pT_K_step0");
991 h50->Write("cT_step0");
992 h60->Write("dca_step0");
993 h70->Write("d0_pi_step0");
994 h80->Write("d0_K_step0");
995 h90->Write("d0xd0_step0");
996 h100->Write("cosPointingAngle_step0");
e8158e56 997 h110->Write("phi_step0");
844d235c 998
999 h01->Write("pT_D0_step1");
1000 h11->Write("rapidity_step1");
1001 h21->Write("cosThetaStar_step1");
1002 h31->Write("pT_pi_step1");
1003 h41->Write("pT_K_step1");
1004 h51->Write("cT_step1");
1005 h61->Write("dca_step1");
1006 h71->Write("d0_pi_step1");
1007 h81->Write("d0_K_step1");
1008 h91->Write("d0xd0_step1");
1009 h101->Write("cosPointingAngle_step1");
e8158e56 1010 h111->Write("phi_step1");
844d235c 1011
1012 h02->Write("pT_D0_step2");
1013 h12->Write("rapidity_step2");
1014 h22->Write("cosThetaStar_step2");
1015 h32->Write("pT_pi_step2");
1016 h42->Write("pT_K_step2");
1017 h52->Write("cT_step2");
1018 h62->Write("dca_step2");
1019 h72->Write("d0_pi_step2");
1020 h80->Write("d0_K_step2");
1021 h92->Write("d0xd0_step2");
1022 h102->Write("cosPointingAngle_step2");
e8158e56 1023 h112->Write("phi_step2");
844d235c 1024
1025 file_projection->Close();
1026
a5d92cb4 1027
1028
f5ec6c30 1029 /*
1030 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1031 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1032 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1033 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1034
1035 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1036 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1037 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1038 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1039 */
1040}
1041
1042//___________________________________________________________________________
1043void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1044 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1045 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1046 //
1047 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1048
1049 //slot #1
1050 OpenFile(1);
1051 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
1052}
1053
1054//___________________________________________________________________________
1055Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1056
1057 //
1058 // to calculate cos(ThetaStar) for generated particle
1059 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1060 // (see where the function is called)
1061 //
1062
1063 Int_t pdgvtx = mcPart->GetPdgCode();
1064 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1065 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1066 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1067 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1068 AliDebug(2,"This is a D0");
1069 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1070 mcPartDaughter0 = mcPartDaughter1;
1071 mcPartDaughter1 = mcPartdummy;
1072 }
1073 else{
1074 AliInfo("D0bar");
1075 }
1076 */
1077 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1078 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1079 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1080 AliDebug(2,"D0");
1081 }
1082 else{
1083 AliDebug(2,"D0bar");
1084 }
1085 if (pdgprong0 == 211){
1086 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1087 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1088 mcPartDaughter0 = mcPartDaughter1;
1089 mcPartDaughter1 = mcPartdummy;
1090 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1091 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1092 }
1093
1094 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1095 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1096 Double_t massp[2];
1097 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1098 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1099
1100 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);
1101
1102 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1103 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1104 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1105 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1106 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1107 Double_t beta = p/e;
1108 Double_t gamma = e/massvtx;
1109 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1110 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1111
1112 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1113
1114 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1115 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1116 AliDebug(2,Form("cts = %f", cts));
1117 return cts;
1118}
1119//___________________________________________________________________________
1120Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1121
1122 //
1123 // to calculate cT for generated particle
1124 //
1125
1126 Double_t xmcPart[3] = {0,0,0};
1127 Double_t xdaughter0[3] = {0,0,0};
1128 Double_t xdaughter1[3] = {0,0,0};
1129 mcPart->XvYvZv(xmcPart); // cm
1130 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1131 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1132 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1133 xmcPart[1]*xmcPart[1]+
1134 xmcPart[2]*xmcPart[2]);
1135 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1136 xdaughter0[1]*xdaughter0[1]+
1137 xdaughter0[2]*xdaughter0[2]);
1138 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1139 xdaughter1[1]*xdaughter1[1]+
1140 xdaughter1[2]*xdaughter1[2]);
1141
1142 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));
1143 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));
1144 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));
1145
1146 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1147 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1148 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1149
1150 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1151 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1152 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1153
e8158e56 1154 if ((cT0 - cT1)>1E-5) {
1155 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 1156 }
1157
1158 // calculating cT from cT0
1159
1160 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1161 Double_t p = mcPart-> P();
1162 Double_t cT = cT0*mass/p;
1163 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1164 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1165 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1166 return cT;
1167}
1168//________________________________________________________________________________
1169Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1170
1171 //
1172 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1173 //
1174
1175 Bool_t isOk = kFALSE;
1176
1177 // check whether the D0 decays in pi+K
1178 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 // could use a cut in the CF, but not clear how to define a TDecayChannel
1180 // implemented for the time being as a cut on the number of daughters - see later when
1181 // getting the daughters
1182
1183 // getting the daughters
1184 Int_t daughter0 = mcPart->GetDaughter(0);
1185 Int_t daughter1 = mcPart->GetDaughter(1);
1186 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1187 if (daughter0 == 0 || daughter1 == 0) {
1188 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1189 return isOk;
1190 }
421623bd 1191 if (TMath::Abs(daughter1 - daughter0) != 1) {
f5ec6c30 1192 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1193 return isOk;
1194 }
1195 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1196 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1197 if (!mcPartDaughter0 || !mcPartDaughter1) {
1198 AliWarning("At least one Daughter Particle not found in tree, skipping");
1199 return isOk;
1200 }
2bcc508a 1201 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1202 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1203 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1204 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1205 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1206 return isOk;
1207 }
1208
f5ec6c30 1209 Double_t vtx1[3] = {0,0,0}; // primary vertex
1210 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1211 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1212 mcPart->XvYvZv(vtx1); // cm
1213 // getting vertex from daughters
1214 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1215 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1216 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1217 AliError("Daughters have different secondary vertex, skipping the track");
1218 return isOk;
1219 }
1220 Int_t nprongs = 2;
1221 Short_t charge = 0;
1222 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1223 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1224 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1225 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1226 // inverting in case the positive daughter is the second one
1227 positiveDaugh = mcPartDaughter1;
1228 negativeDaugh = mcPartDaughter0;
1229 }
1230 // getting the momentum from the daughters
1231 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1232 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1233 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1234
1235 Double_t d0[2] = {0.,0.};
1236
1237 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1238
1239 Double_t cosThetaStar = 0.;
1240 Double_t cosThetaStarD0 = 0.;
1241 Double_t cosThetaStarD0bar = 0.;
1242 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1243 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1244 if (mcPart->GetPdgCode() == 421){ // D0
1245 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1246 cosThetaStar = cosThetaStarD0;
1247 }
1248 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1249 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1250 cosThetaStar = cosThetaStarD0bar;
1251 }
1252 else{
1253 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1254 return vectorMC;
1255 }
1256 if (cosThetaStar < -1 || cosThetaStar > 1) {
1257 AliWarning("Invalid value for cosine Theta star, returning");
1258 return isOk;
1259 }
1260
1261 // calculate cos(Theta*) according to the method implemented herein
1262
1263 Double_t cts = 9999.;
1264 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1265 if (cts < -1 || cts > 1) {
1266 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1267 }
1268 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1269 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1270 }
1271
1272 Double_t cT = decay->Ct(421);
1273
1274 // calculate cT from the method implemented herein
1275 Double_t cT1 = 0.;
1276 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1277
1278 if (TMath::Abs(cT1 - cT)>0.001) {
1279 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1280 }
1281
1282 // get the pT of the daughters
1283
1284 Double_t pTpi = 0.;
1285 Double_t pTK = 0.;
1286
1287 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1288 pTpi = mcPartDaughter0->Pt();
1289 pTK = mcPartDaughter1->Pt();
1290 }
1291 else {
1292 pTpi = mcPartDaughter1->Pt();
1293 pTK = mcPartDaughter0->Pt();
1294 }
1295
1296 vectorMC[0] = mcPart->Pt();
1297 vectorMC[1] = mcPart->Y() ;
1298 vectorMC[2] = cosThetaStar ;
1299 vectorMC[3] = pTpi ;
1300 vectorMC[4] = pTK ;
1301 vectorMC[5] = cT*1.E4 ; // in micron
e8158e56 1302 vectorMC[6] = mcPart->Phi() ;
f5ec6c30 1303 isOk = kTRUE;
1304 return isOk;
1305}
e8158e56 1306//_________________________________________________________________________________________________
1307Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1308
1309 //
1310 // checking whether the very mother of the D0 is a charm or a bottom quark
1311 //
1312
1313 Int_t pdgGranma = 0;
1314 Int_t mother = 0;
1315 mother = mcPart->GetMother();
1316 Int_t istep = 0;
1317 while (mother >0 ){
1318 istep++;
1319 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1320 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1321 pdgGranma = mcGranma->GetPdgCode();
1322 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1323 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1324 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1325 break;
1326 }
1327 mother = mcGranma->GetMother();
1328 }
1329 return pdgGranma;
1330}
f5ec6c30 1331