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