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