1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------------
17 // Class for HF corrections as a function of many variables
18 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19 // Reco Acc + ITS Cl + PPR cuts
20 // 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
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 //-----------------------------------------------------------------------
32 #include <TParticle.h>
33 #include <TDatabasePDG.h>
38 #include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
40 #include "AliMCEvent.h"
41 #include "AliCFManager.h"
42 #include "AliCFContainer.h"
44 #include "AliAODEvent.h"
45 #include "AliAODRecoDecay.h"
46 #include "AliAODRecoDecayHF.h"
47 #include "AliAODRecoDecayHF2Prong.h"
48 #include "AliAODMCParticle.h"
49 #include "AliESDtrack.h"
51 #include "THnSparse.h"
53 //__________________________________________________________________________
54 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
58 fHistEventsProcessed(0x0),
66 fCountRecoITSClusters(0),
69 fFillFromGenerated(kFALSE),
77 //___________________________________________________________________________
78 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
79 AliAnalysisTaskSE(name),
82 fHistEventsProcessed(0x0),
90 fCountRecoITSClusters(0),
93 fFillFromGenerated(kFALSE),
98 // Constructor. Initialization of Inputs and Outputs
100 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
102 DefineInput(0) and DefineOutput(0)
103 are taken care of by AliAnalysisTaskSE constructor
105 DefineOutput(1,TH1I::Class());
106 DefineOutput(2,AliCFContainer::Class());
107 DefineOutput(3,THnSparseD::Class());
110 //___________________________________________________________________________
111 AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
114 // Assignment operator
117 AliAnalysisTaskSE::operator=(c) ;
119 fCFManager = c.fCFManager;
120 fHistEventsProcessed = c.fHistEventsProcessed;
125 //___________________________________________________________________________
126 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
127 AliAnalysisTaskSE(c),
129 fCFManager(c.fCFManager),
130 fHistEventsProcessed(c.fHistEventsProcessed),
131 fCorrelation(c.fCorrelation),
132 fCountMC(c.fCountMC),
133 fCountAcc(c.fCountAcc),
134 fCountVertex(c.fCountVertex),
135 fCountRefit(c.fCountRefit),
136 fCountReco(c.fCountReco),
137 fCountRecoAcc(c.fCountRecoAcc),
138 fCountRecoITSClusters(c.fCountRecoITSClusters),
139 fCountRecoPPR(c.fCountRecoPPR),
141 fFillFromGenerated(c.fFillFromGenerated),
142 fMinITSClusters(c.fMinITSClusters),
143 fAcceptanceUnf(c.fAcceptanceUnf)
150 //___________________________________________________________________________
151 AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
155 if (fCFManager) delete fCFManager ;
156 if (fHistEventsProcessed) delete fHistEventsProcessed ;
157 if (fCorrelation) delete fCorrelation ;
160 //_________________________________________________
161 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
164 // Main loop function
167 if (fFillFromGenerated){
168 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
172 Error("UserExec","NO EVENT FOUND!");
177 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
178 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
179 fCFManager->SetEventInfo(aodEvent);
181 // MC-event selection
182 Double_t containerInput[12] ;
183 Double_t containerInputMC[12] ;
185 //loop on the MC event
187 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
188 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
191 Int_t icountReco = 0;
192 Int_t icountVertex = 0;
193 Int_t icountRefit = 0;
194 Int_t icountRecoAcc = 0;
195 Int_t icountRecoITSClusters = 0;
196 Int_t icountRecoPPR = 0;
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;
206 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
207 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
208 if (mcPart->GetPdgCode() == 4) cquarks++;
209 if (mcPart->GetPdgCode() == -4) cquarks++;
211 AliWarning("Particle not found in tree, skipping");
215 // check the MC-level cuts
217 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
218 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
219 Int_t abspdgGranma = TMath::Abs(pdgGranma);
220 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
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));
222 continue; // skipping particles that don't come from c quark
225 // if (TMath::Abs(pdgGranma)!=4) {
227 // fill the container for Gen-level selection
228 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
229 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
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);
245 // check the MC-Acceptance level cuts
246 // since standard CF functions are not applicable, using Kine Cuts on daughters
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...");
254 if (TMath::Abs(daughter1 - daughter0) != 1) {
255 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
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...");
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)) {
276 AliInfo("Inconsistency with CF cut for daughter 0!");
278 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
279 AliInfo("Inconsistency with CF cut for daughter 1!");
281 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
284 // check on the vertex
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");
291 // check on the kTPCrefit and kITSrefit conditions of the daughters
292 Bool_t refitFlag = kFALSE;
293 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
294 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
295 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1) && (((track->GetStatus()&AliESDtrack::kTPCrefit)==AliESDtrack::kTPCrefit) && ((track->GetStatus()&AliESDtrack::kITSrefit)==AliESDtrack::kITSrefit))){
300 printf("Refit cut passed\n");
301 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
305 AliDebug(3,"Refit cut not passed\n");
309 AliDebug(3,"Vertex cut not passed\n");
313 AliDebug(3,"Acceptance cut not passed\n");
317 AliDebug(3,"Problems in filling the container");
322 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
324 AliInfo(Form("Found %i MC particles that are D0!!",icountMC));
325 AliInfo(Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
326 AliInfo(Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
327 AliInfo(Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
329 // Now go to rec level
330 fCountMC += icountMC;
331 fCountAcc += icountAcc;
333 // load heavy flavour vertices
334 TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
335 if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
336 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
338 Int_t pdgDgD0toKpi[2]={321,211};
340 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
342 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
343 Bool_t unsetvtx=kFALSE;
344 if(!d0tokpi->GetOwnPrimaryVtx()) {
345 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
349 // find associated MC particle
350 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
353 AliDebug(2,"No MC particle found");
356 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
358 AliWarning("Could not find associated MC in AOD MC tree");
361 // check whether the daughters have kTPCrefit and kITSrefit set
362 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
363 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
364 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
365 // skipping if at least one daughter does not have kTPCrefit or kITSrefit
369 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
371 // skipping cause vertex is not ok
375 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
377 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
378 // skipping candidate
382 // check if associated MC v0 passes the cuts
383 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
384 AliDebug(2, "Skipping the particles due to cuts");
387 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
388 Int_t abspdgGranma = TMath::Abs(pdgGranma);
389 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
390 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));
391 continue; // skipping particles that don't come from c quark
394 // fill the container...
395 Double_t pt = d0tokpi->Pt();
396 Double_t rapidity = d0tokpi->YD0();
398 Double_t cosThetaStar = 9999.;
401 Double_t dca = d0tokpi->GetDCA();
404 Double_t d0xd0 = d0tokpi->Prodd0d0();
405 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
406 Double_t phi = d0tokpi->Phi();
407 Int_t pdgCode = mcVtxHF->GetPdgCode();
409 cosThetaStar = d0tokpi->CosThetaStarD0();
410 pTpi = d0tokpi->PtProng(0);
411 pTK = d0tokpi->PtProng(1);
412 d0pi = d0tokpi->Getd0Prong(0);
413 d0K = d0tokpi->Getd0Prong(1);
416 cosThetaStar = d0tokpi->CosThetaStarD0bar();
417 pTpi = d0tokpi->PtProng(1);
418 pTK = d0tokpi->PtProng(0);
419 d0pi = d0tokpi->Getd0Prong(1);
420 d0K = d0tokpi->Getd0Prong(0);
423 Double_t cT = d0tokpi->CtD0();
425 if (!fFillFromGenerated){
426 // ...either with reconstructed values....
427 containerInput[0] = pt;
428 containerInput[1] = rapidity;
429 containerInput[2] = cosThetaStar;
430 containerInput[3] = pTpi;
431 containerInput[4] = pTK;
432 containerInput[5] = cT*1.E4; // in micron
433 containerInput[6] = dca*1.E4; // in micron
434 containerInput[7] = d0pi*1.E4; // in micron
435 containerInput[8] = d0K*1.E4; // in micron
436 containerInput[9] = d0xd0*1.E8; // in micron^2
437 containerInput[10] = cosPointingAngle; // in micron
438 containerInput[11] = phi;
441 // ... or with generated values
442 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
443 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
444 containerInput[0] = vectorMC[0];
445 containerInput[1] = vectorMC[1] ;
446 containerInput[2] = vectorMC[2] ;
447 containerInput[3] = vectorMC[3] ;
448 containerInput[4] = vectorMC[4] ;
449 containerInput[5] = vectorMC[5] ; // in micron
450 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
451 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
452 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
453 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
454 containerInput[10] = 1.01; // dummy value, meaningless in MC
455 containerInput[11] = vectorMC[6];
458 AliDebug(3,"Problems in filling the container");
462 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]));
464 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
467 Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
468 Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
469 if (acceptanceProng0 && acceptanceProng1) {
470 AliDebug(2,"D0 reco daughters in acceptance");
471 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
476 Double_t fill[4]; //fill response matrix
478 // dimensions 0&1 : pt,eta (Rec)
483 // dimensions 2&3 : pt,eta (MC)
485 fill[2] = mcVtxHF->Pt();
486 fill[3] = mcVtxHF->Y();
488 fCorrelation->Fill(fill);
492 // cut on the min n. of clusters in ITS
494 for(Int_t l=0;l<6;l++) if(TESTBIT(d0tokpi->GetITSClusterMap(),l)) ncls0++;
495 AliDebug(2, Form("n clusters = %d", ncls0));
496 if (ncls0 >= fMinITSClusters){
497 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
498 icountRecoITSClusters++;
499 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));
502 Double_t cuts[6] = {9999999., 1.1, 0., 9999999., 9999999., 0.};
511 else if (pt > 1 && pt <= 2){
519 else if (pt > 2 && pt <= 3){
527 else if (pt > 3 && pt <= 5){
543 if (dca*1E4 < cuts[0]
544 && TMath::Abs(cosThetaStar) < cuts[1]
547 && TMath::Abs(d0pi*1E4) < cuts[3]
548 && TMath::Abs(d0K*1E4) < cuts[3]
549 && d0xd0*1E8 < cuts[4]
550 && cosPointingAngle > cuts[5]
553 AliDebug(2,"Particle passed PPR cuts");
554 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
557 if(!fAcceptanceUnf){ // unfolding
559 Double_t fill[4]; //fill response matrix
561 // dimensions 0&1 : pt,eta (Rec)
566 // dimensions 2&3 : pt,eta (MC)
568 fill[2] = mcVtxHF->Pt();
569 fill[3] = mcVtxHF->Y();
571 fCorrelation->Fill(fill);
576 AliDebug(2,"Particle skipped due to PPR cuts");
577 if (dca*1E4 > cuts[0]){
578 AliDebug(2,"Problems with dca");
580 if ( TMath::Abs(cosThetaStar) > cuts[1]){
581 AliDebug(2,"Problems with cosThetaStar");
584 AliDebug(2,"Problems with pTpi");
587 AliDebug(2,"Problems with pTK");
589 if (TMath::Abs(d0pi*1E4) > cuts[3]){
590 AliDebug(2,"Problems with d0pi");
592 if (TMath::Abs(d0K*1E4) > cuts[3]){
593 AliDebug(2,"Problems with d0K");
595 if (d0xd0*1E8 > cuts[4]){
596 AliDebug(2,"Problems with d0xd0");
598 if (cosPointingAngle < cuts[5]){
599 AliDebug(2,"Problems with cosPointingAngle");
605 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
606 } // end loop on D0->Kpi
608 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
610 fCountReco+= icountReco;
611 fCountVertex+= icountVertex;
612 fCountRefit+= icountRefit;
613 fCountRecoAcc+= icountRecoAcc;
614 fCountRecoITSClusters+= icountRecoITSClusters;
615 fCountRecoPPR+= icountRecoPPR;
617 fHistEventsProcessed->Fill(0);
618 /* PostData(0) is taken care of by AliAnalysisTaskSE */
619 PostData(1,fHistEventsProcessed) ;
620 PostData(2,fCFManager->GetParticleContainer()) ;
621 PostData(3,fCorrelation) ;
625 //___________________________________________________________________________
626 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
628 // The Terminate() function is the last function to be called during
629 // a query. It always runs on the client, it can be used to present
630 // the results graphically or save the results to file.
632 AliAnalysisTaskSE::Terminate();
634 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
635 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
636 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
637 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));
638 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
639 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));
640 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));
641 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
643 // draw some example plots....
645 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
647 printf("CONTAINER NOT FOUND\n");
650 // projecting the containers to obtain histograms
651 // first argument = variable, second argument = step
654 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
655 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
656 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
657 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
658 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
659 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
660 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
661 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
662 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
663 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
664 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
665 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
667 // MC-Acceptance level
668 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
669 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
670 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
671 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
672 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
673 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
674 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
675 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
676 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
677 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
678 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
679 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
682 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
683 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
684 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
685 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
686 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
687 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
688 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
689 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
690 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
691 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
692 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
693 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
695 h00->SetTitle("pT_D0 (GeV/c)");
696 h10->SetTitle("rapidity");
697 h20->SetTitle("cosThetaStar");
698 h30->SetTitle("pT_pi (GeV/c)");
699 h40->SetTitle("pT_K (Gev/c)");
700 h50->SetTitle("cT (#mum)");
701 h60->SetTitle("dca (#mum)");
702 h70->SetTitle("d0_pi (#mum)");
703 h80->SetTitle("d0_K (#mum)");
704 h90->SetTitle("d0xd0 (#mum^2)");
705 h100->SetTitle("cosPointingAngle");
706 h100->SetTitle("phi (rad)");
708 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
709 h10->GetXaxis()->SetTitle("rapidity");
710 h20->GetXaxis()->SetTitle("cosThetaStar");
711 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
712 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
713 h50->GetXaxis()->SetTitle("cT (#mum)");
714 h60->GetXaxis()->SetTitle("dca (#mum)");
715 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
716 h80->GetXaxis()->SetTitle("d0_K (#mum)");
717 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
718 h100->GetXaxis()->SetTitle("cosPointingAngle");
719 h110->GetXaxis()->SetTitle("phi (rad)");
721 h01->SetTitle("pT_D0 (GeV/c)");
722 h11->SetTitle("rapidity");
723 h21->SetTitle("cosThetaStar");
724 h31->SetTitle("pT_pi (GeV/c)");
725 h41->SetTitle("pT_K (Gev/c)");
726 h51->SetTitle("cT (#mum)");
727 h61->SetTitle("dca (#mum)");
728 h71->SetTitle("d0_pi (#mum)");
729 h81->SetTitle("d0_K (#mum)");
730 h91->SetTitle("d0xd0 (#mum^2)");
731 h101->SetTitle("cosPointingAngle");
732 h111->GetXaxis()->SetTitle("phi (rad)");
734 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
735 h11->GetXaxis()->SetTitle("rapidity");
736 h21->GetXaxis()->SetTitle("cosThetaStar");
737 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
738 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
739 h51->GetXaxis()->SetTitle("cT (#mum)");
740 h61->GetXaxis()->SetTitle("dca (#mum)");
741 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
742 h81->GetXaxis()->SetTitle("d0_K (#mum)");
743 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
744 h101->GetXaxis()->SetTitle("cosPointingAngle");
745 h111->GetXaxis()->SetTitle("phi (rad)");
747 h04->SetTitle("pT_D0 (GeV/c)");
748 h14->SetTitle("rapidity");
749 h24->SetTitle("cosThetaStar");
750 h34->SetTitle("pT_pi (GeV/c)");
751 h44->SetTitle("pT_K (Gev/c)");
752 h54->SetTitle("cT (#mum)");
753 h64->SetTitle("dca (#mum)");
754 h74->SetTitle("d0_pi (#mum)");
755 h84->SetTitle("d0_K (#mum)");
756 h94->SetTitle("d0xd0 (#mum^2)");
757 h104->SetTitle("cosPointingAngle");
758 h114->GetXaxis()->SetTitle("phi (rad)");
760 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
761 h14->GetXaxis()->SetTitle("rapidity");
762 h24->GetXaxis()->SetTitle("cosThetaStar");
763 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
764 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
765 h54->GetXaxis()->SetTitle("cT (#mum)");
766 h64->GetXaxis()->SetTitle("dca (#mum)");
767 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
768 h84->GetXaxis()->SetTitle("d0_K (#mum)");
769 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
770 h104->GetXaxis()->SetTitle("cosPointingAngle");
771 h114->GetXaxis()->SetTitle("phi (rad)");
773 Double_t max0 = h00->GetMaximum();
774 Double_t max1 = h10->GetMaximum();
775 Double_t max2 = h20->GetMaximum();
776 Double_t max3 = h30->GetMaximum();
777 Double_t max4 = h40->GetMaximum();
778 Double_t max5 = h50->GetMaximum();
779 Double_t max6 = h60->GetMaximum();
780 Double_t max7 = h70->GetMaximum();
781 Double_t max8 = h80->GetMaximum();
782 Double_t max9 = h90->GetMaximum();
783 Double_t max10 = h100->GetMaximum();
784 Double_t max11 = h110->GetMaximum();
786 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
787 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
788 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
789 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
790 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
791 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
792 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
793 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
794 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
795 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
796 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
797 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
799 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
800 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
801 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
802 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
803 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
804 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
805 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
806 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
807 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
808 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
809 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
810 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
812 h00->SetMarkerStyle(20);
813 h10->SetMarkerStyle(24);
814 h20->SetMarkerStyle(21);
815 h30->SetMarkerStyle(25);
816 h40->SetMarkerStyle(27);
817 h50->SetMarkerStyle(28);
818 h60->SetMarkerStyle(20);
819 h70->SetMarkerStyle(24);
820 h80->SetMarkerStyle(21);
821 h90->SetMarkerStyle(25);
822 h100->SetMarkerStyle(27);
823 h110->SetMarkerStyle(28);
825 h00->SetMarkerColor(2);
826 h10->SetMarkerColor(2);
827 h20->SetMarkerColor(2);
828 h30->SetMarkerColor(2);
829 h40->SetMarkerColor(2);
830 h50->SetMarkerColor(2);
831 h60->SetMarkerColor(2);
832 h70->SetMarkerColor(2);
833 h80->SetMarkerColor(2);
834 h90->SetMarkerColor(2);
835 h100->SetMarkerColor(2);
836 h110->SetMarkerColor(2);
838 h01->SetMarkerStyle(20) ;
839 h11->SetMarkerStyle(24) ;
840 h21->SetMarkerStyle(21) ;
841 h31->SetMarkerStyle(25) ;
842 h41->SetMarkerStyle(27) ;
843 h51->SetMarkerStyle(28) ;
844 h61->SetMarkerStyle(20);
845 h71->SetMarkerStyle(24);
846 h81->SetMarkerStyle(21);
847 h91->SetMarkerStyle(25);
848 h101->SetMarkerStyle(27);
849 h111->SetMarkerStyle(28);
851 h01->SetMarkerColor(8);
852 h11->SetMarkerColor(8);
853 h21->SetMarkerColor(8);
854 h31->SetMarkerColor(8);
855 h41->SetMarkerColor(8);
856 h51->SetMarkerColor(8);
857 h61->SetMarkerColor(8);
858 h71->SetMarkerColor(8);
859 h81->SetMarkerColor(8);
860 h91->SetMarkerColor(8);
861 h101->SetMarkerColor(8);
862 h111->SetMarkerColor(8);
864 h04->SetMarkerStyle(20) ;
865 h14->SetMarkerStyle(24) ;
866 h24->SetMarkerStyle(21) ;
867 h34->SetMarkerStyle(25) ;
868 h44->SetMarkerStyle(27) ;
869 h54->SetMarkerStyle(28) ;
870 h64->SetMarkerStyle(20);
871 h74->SetMarkerStyle(24);
872 h84->SetMarkerStyle(21);
873 h94->SetMarkerStyle(25);
874 h104->SetMarkerStyle(27);
875 h114->SetMarkerStyle(28);
877 h04->SetMarkerColor(4);
878 h14->SetMarkerColor(4);
879 h24->SetMarkerColor(4);
880 h34->SetMarkerColor(4);
881 h44->SetMarkerColor(4);
882 h54->SetMarkerColor(4);
883 h64->SetMarkerColor(4);
884 h74->SetMarkerColor(4);
885 h84->SetMarkerColor(4);
886 h94->SetMarkerColor(4);
887 h104->SetMarkerColor(4);
888 h114->SetMarkerColor(4);
890 gStyle->SetCanvasColor(0);
891 gStyle->SetFrameFillColor(0);
892 gStyle->SetTitleFillColor(0);
893 gStyle->SetStatColor(0);
895 // drawing in 2 separate canvas for a matter of clearity
896 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
928 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
959 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
990 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1021 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1023 TH2D* corr1 =hcorr->Projection(0,2);
1024 TH2D* corr2 = hcorr->Projection(1,3);
1026 TCanvas * c7 =new TCanvas("c7","",800,400);
1029 corr1->Draw("text");
1031 corr2->Draw("text");
1034 TFile* file_projection = new TFile("file_projection.root","RECREATE");
1038 h00->Write("pT_D0_step0");
1039 h10->Write("rapidity_step0");
1040 h20->Write("cosThetaStar_step0");
1041 h30->Write("pT_pi_step0");
1042 h40->Write("pT_K_step0");
1043 h50->Write("cT_step0");
1044 h60->Write("dca_step0");
1045 h70->Write("d0_pi_step0");
1046 h80->Write("d0_K_step0");
1047 h90->Write("d0xd0_step0");
1048 h100->Write("cosPointingAngle_step0");
1049 h110->Write("phi_step0");
1051 h01->Write("pT_D0_step1");
1052 h11->Write("rapidity_step1");
1053 h21->Write("cosThetaStar_step1");
1054 h31->Write("pT_pi_step1");
1055 h41->Write("pT_K_step1");
1056 h51->Write("cT_step1");
1057 h61->Write("dca_step1");
1058 h71->Write("d0_pi_step1");
1059 h81->Write("d0_K_step1");
1060 h91->Write("d0xd0_step1");
1061 h101->Write("cosPointingAngle_step1");
1062 h111->Write("phi_step1");
1064 h04->Write("pT_D0_step2");
1065 h14->Write("rapidity_step2");
1066 h24->Write("cosThetaStar_step2");
1067 h34->Write("pT_pi_step2");
1068 h44->Write("pT_K_step2");
1069 h54->Write("cT_step2");
1070 h64->Write("dca_step2");
1071 h74->Write("d0_pi_step2");
1072 h80->Write("d0_K_step2");
1073 h94->Write("d0xd0_step2");
1074 h104->Write("cosPointingAngle_step2");
1075 h114->Write("phi_step2");
1077 file_projection->Close();
1082 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1083 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1084 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1085 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1087 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1088 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1089 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1090 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1094 //___________________________________________________________________________
1095 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1096 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1097 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1099 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1103 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
1106 //___________________________________________________________________________
1107 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1110 // to calculate cos(ThetaStar) for generated particle
1111 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1112 // (see where the function is called)
1115 Int_t pdgvtx = mcPart->GetPdgCode();
1116 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1117 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1118 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1119 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1120 AliDebug(2,"This is a D0");
1121 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1122 mcPartDaughter0 = mcPartDaughter1;
1123 mcPartDaughter1 = mcPartdummy;
1129 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1130 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1131 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1135 AliDebug(2,"D0bar");
1137 if (pdgprong0 == 211){
1138 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1139 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1140 mcPartDaughter0 = mcPartDaughter1;
1141 mcPartDaughter1 = mcPartdummy;
1142 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1143 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1146 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1147 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1149 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1150 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1152 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);
1154 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1155 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1156 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1157 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1158 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1159 Double_t beta = p/e;
1160 Double_t gamma = e/massvtx;
1161 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1162 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1164 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1166 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1167 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1168 AliDebug(2,Form("cts = %f", cts));
1171 //___________________________________________________________________________
1172 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1175 // to calculate cT for generated particle
1178 Double_t xmcPart[3] = {0,0,0};
1179 Double_t xdaughter0[3] = {0,0,0};
1180 Double_t xdaughter1[3] = {0,0,0};
1181 mcPart->XvYvZv(xmcPart); // cm
1182 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1183 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1184 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1185 xmcPart[1]*xmcPart[1]+
1186 xmcPart[2]*xmcPart[2]);
1187 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1188 xdaughter0[1]*xdaughter0[1]+
1189 xdaughter0[2]*xdaughter0[2]);
1190 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1191 xdaughter1[1]*xdaughter1[1]+
1192 xdaughter1[2]*xdaughter1[2]);
1194 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));
1195 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));
1196 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));
1198 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1199 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1200 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1202 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1203 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1204 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1206 if ((cT0 - cT1)>1E-5) {
1207 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1210 // calculating cT from cT0
1212 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1213 Double_t p = mcPart-> P();
1214 Double_t cT = cT0*mass/p;
1215 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1216 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1217 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1220 //________________________________________________________________________________
1221 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1224 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1227 Bool_t isOk = kFALSE;
1229 // check whether the D0 decays in pi+K
1230 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1231 // could use a cut in the CF, but not clear how to define a TDecayChannel
1232 // implemented for the time being as a cut on the number of daughters - see later when
1233 // getting the daughters
1235 // getting the daughters
1236 Int_t daughter0 = mcPart->GetDaughter(0);
1237 Int_t daughter1 = mcPart->GetDaughter(1);
1238 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1239 if (daughter0 == 0 || daughter1 == 0) {
1240 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1243 if (TMath::Abs(daughter1 - daughter0) != 1) {
1244 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1247 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1248 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1249 if (!mcPartDaughter0 || !mcPartDaughter1) {
1250 AliWarning("At least one Daughter Particle not found in tree, skipping");
1253 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1254 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1255 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1256 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1257 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1261 Double_t vtx1[3] = {0,0,0}; // primary vertex
1262 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1263 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1264 mcPart->XvYvZv(vtx1); // cm
1265 // getting vertex from daughters
1266 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1267 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1268 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1269 AliError("Daughters have different secondary vertex, skipping the track");
1274 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1275 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1276 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1277 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1278 // inverting in case the positive daughter is the second one
1279 positiveDaugh = mcPartDaughter1;
1280 negativeDaugh = mcPartDaughter0;
1282 // getting the momentum from the daughters
1283 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1284 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1285 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1287 Double_t d0[2] = {0.,0.};
1289 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1291 Double_t cosThetaStar = 0.;
1292 Double_t cosThetaStarD0 = 0.;
1293 Double_t cosThetaStarD0bar = 0.;
1294 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1295 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1296 if (mcPart->GetPdgCode() == 421){ // D0
1297 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1298 cosThetaStar = cosThetaStarD0;
1300 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1301 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1302 cosThetaStar = cosThetaStarD0bar;
1305 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1308 if (cosThetaStar < -1 || cosThetaStar > 1) {
1309 AliWarning("Invalid value for cosine Theta star, returning");
1313 // calculate cos(Theta*) according to the method implemented herein
1315 Double_t cts = 9999.;
1316 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1317 if (cts < -1 || cts > 1) {
1318 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1320 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1321 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1324 Double_t cT = decay->Ct(421);
1326 // calculate cT from the method implemented herein
1328 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1330 if (TMath::Abs(cT1 - cT)>0.001) {
1331 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1334 // get the pT of the daughters
1339 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1340 pTpi = mcPartDaughter0->Pt();
1341 pTK = mcPartDaughter1->Pt();
1344 pTpi = mcPartDaughter1->Pt();
1345 pTK = mcPartDaughter0->Pt();
1348 vectorMC[0] = mcPart->Pt();
1349 vectorMC[1] = mcPart->Y() ;
1350 vectorMC[2] = cosThetaStar ;
1351 vectorMC[3] = pTpi ;
1353 vectorMC[5] = cT*1.E4 ; // in micron
1354 vectorMC[6] = mcPart->Phi() ;
1358 //_________________________________________________________________________________________________
1359 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1362 // checking whether the very mother of the D0 is a charm or a bottom quark
1365 Int_t pdgGranma = 0;
1367 mother = mcPart->GetMother();
1371 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1372 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1373 pdgGranma = mcGranma->GetPdgCode();
1374 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1375 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1376 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1379 mother = mcGranma->GetMother();