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) AliDebug(2,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 AliDebug(2,"Inconsistency with CF cut for daughter 0!");
278 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
279 AliDebug(2,"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 = kTRUE;
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)) {
296 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
302 printf("Refit cut passed\n");
303 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
307 AliDebug(3,"Refit cut not passed\n");
311 AliDebug(3,"Vertex cut not passed\n");
315 AliDebug(3,"Acceptance cut not passed\n");
319 AliDebug(3,"Problems in filling the container");
324 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
326 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
327 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
328 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
329 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
331 // Now go to rec level
332 fCountMC += icountMC;
333 fCountAcc += icountAcc;
335 // load heavy flavour vertices
336 TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
337 if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
338 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
340 Int_t pdgDgD0toKpi[2]={321,211};
342 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
344 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
345 Bool_t unsetvtx=kFALSE;
346 if(!d0tokpi->GetOwnPrimaryVtx()) {
347 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
351 // find associated MC particle
352 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
355 AliDebug(2,"No MC particle found");
358 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
360 AliWarning("Could not find associated MC in AOD MC tree");
363 // check whether the daughters have kTPCrefit and kITSrefit set
364 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
365 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
366 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
367 // skipping if at least one daughter does not have kTPCrefit or kITSrefit
371 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
373 // skipping cause vertex is not ok
377 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
379 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
380 // skipping candidate
384 // check if associated MC v0 passes the cuts
385 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
386 AliDebug(2, "Skipping the particles due to cuts");
389 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
390 Int_t abspdgGranma = TMath::Abs(pdgGranma);
391 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
392 AliDebug(2,Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
393 continue; // skipping particles that don't come from c quark
396 // fill the container...
397 Double_t pt = d0tokpi->Pt();
398 Double_t rapidity = d0tokpi->YD0();
400 Double_t cosThetaStar = 9999.;
403 Double_t dca = d0tokpi->GetDCA();
406 Double_t d0xd0 = d0tokpi->Prodd0d0();
407 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
408 Double_t phi = d0tokpi->Phi();
409 Int_t pdgCode = mcVtxHF->GetPdgCode();
411 cosThetaStar = d0tokpi->CosThetaStarD0();
412 pTpi = d0tokpi->PtProng(0);
413 pTK = d0tokpi->PtProng(1);
414 d0pi = d0tokpi->Getd0Prong(0);
415 d0K = d0tokpi->Getd0Prong(1);
416 invMass=d0tokpi->InvMassD0();
419 cosThetaStar = d0tokpi->CosThetaStarD0bar();
420 pTpi = d0tokpi->PtProng(1);
421 pTK = d0tokpi->PtProng(0);
422 d0pi = d0tokpi->Getd0Prong(1);
423 d0K = d0tokpi->Getd0Prong(0);
424 invMass=d0tokpi->InvMassD0bar();
427 Double_t cT = d0tokpi->CtD0();
429 if (!fFillFromGenerated){
430 // ...either with reconstructed values....
431 containerInput[0] = pt;
432 containerInput[1] = rapidity;
433 containerInput[2] = cosThetaStar;
434 containerInput[3] = pTpi;
435 containerInput[4] = pTK;
436 containerInput[5] = cT*1.E4; // in micron
437 containerInput[6] = dca*1.E4; // in micron
438 containerInput[7] = d0pi*1.E4; // in micron
439 containerInput[8] = d0K*1.E4; // in micron
440 containerInput[9] = d0xd0*1.E8; // in micron^2
441 containerInput[10] = cosPointingAngle; // in micron
442 containerInput[11] = phi;
445 // ... or with generated values
446 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
447 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
448 containerInput[0] = vectorMC[0];
449 containerInput[1] = vectorMC[1] ;
450 containerInput[2] = vectorMC[2] ;
451 containerInput[3] = vectorMC[3] ;
452 containerInput[4] = vectorMC[4] ;
453 containerInput[5] = vectorMC[5] ; // in micron
454 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
455 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
456 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
457 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
458 containerInput[10] = 1.01; // dummy value, meaningless in MC
459 containerInput[11] = vectorMC[6];
462 AliDebug(3,"Problems in filling the container");
466 AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
468 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
471 Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
472 Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
473 if (acceptanceProng0 && acceptanceProng1) {
474 AliDebug(2,"D0 reco daughters in acceptance");
475 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
480 Double_t fill[4]; //fill response matrix
482 // dimensions 0&1 : pt,eta (Rec)
487 // dimensions 2&3 : pt,eta (MC)
489 fill[2] = mcVtxHF->Pt();
490 fill[3] = mcVtxHF->Y();
492 fCorrelation->Fill(fill);
496 // cut on the min n. of clusters in ITS
497 Int_t ncls0=0,ncls1=0;
498 for(Int_t l=0;l<6;l++) {
499 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
500 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
502 AliDebug(2, Form("n clusters = %d", ncls0));
503 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters) {
504 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
505 icountRecoITSClusters++;
506 AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
509 //added mass cut for D0 analysis
510 Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
511 //cuts of the D0 analysis (looser) see AliAnalysisTaskSED0Mass.cxx
513 cuts[0] = 400; //dca (um)
514 cuts[1] = 0.8; //cosTstar
515 cuts[2] = 0.5; //pt (GeV/c)
516 cuts[3] = 500; //d0 (um)
517 cuts[4] = -25000; //d0xd0 (um^2)
518 cuts[5] = 0.7; //cosTpointing
520 else if (pt > 1 && pt <= 3){
528 else if (pt > 3 && pt <= 5){
536 // else if (pt > 3 && pt <= 5){
553 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
554 if (TMath::Abs(invMass-mD0PDG) < cuts[6]
556 && TMath::Abs(cosThetaStar) < cuts[1]
559 && TMath::Abs(d0pi*1E4) < cuts[3]
560 && TMath::Abs(d0K*1E4) < cuts[3]
561 && d0xd0*1E8 < cuts[4]
562 && cosPointingAngle > cuts[5]
565 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
566 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
569 if(!fAcceptanceUnf){ // unfolding
571 Double_t fill[4]; //fill response matrix
573 // dimensions 0&1 : pt,eta (Rec)
578 // dimensions 2&3 : pt,eta (MC)
580 fill[2] = mcVtxHF->Pt();
581 fill[3] = mcVtxHF->Y();
583 fCorrelation->Fill(fill);
588 AliDebug(2,"Particle skipped due to PPR cuts");
589 if (dca*1E4 > cuts[0]){
590 AliDebug(2,"Problems with dca");
592 if ( TMath::Abs(cosThetaStar) > cuts[1]){
593 AliDebug(2,"Problems with cosThetaStar");
596 AliDebug(2,"Problems with pTpi");
599 AliDebug(2,"Problems with pTK");
601 if (TMath::Abs(d0pi*1E4) > cuts[3]){
602 AliDebug(2,"Problems with d0pi");
604 if (TMath::Abs(d0K*1E4) > cuts[3]){
605 AliDebug(2,"Problems with d0K");
607 if (d0xd0*1E8 > cuts[4]){
608 AliDebug(2,"Problems with d0xd0");
610 if (cosPointingAngle < cuts[5]){
611 AliDebug(2,"Problems with cosPointingAngle");
617 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
618 } // end loop on D0->Kpi
620 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
622 fCountReco+= icountReco;
623 fCountVertex+= icountVertex;
624 fCountRefit+= icountRefit;
625 fCountRecoAcc+= icountRecoAcc;
626 fCountRecoITSClusters+= icountRecoITSClusters;
627 fCountRecoPPR+= icountRecoPPR;
629 fHistEventsProcessed->Fill(0);
630 /* PostData(0) is taken care of by AliAnalysisTaskSE */
631 PostData(1,fHistEventsProcessed) ;
632 PostData(2,fCFManager->GetParticleContainer()) ;
633 PostData(3,fCorrelation) ;
637 //___________________________________________________________________________
638 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
640 // The Terminate() function is the last function to be called during
641 // a query. It always runs on the client, it can be used to present
642 // the results graphically or save the results to file.
644 AliAnalysisTaskSE::Terminate();
646 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
647 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
648 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
649 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fEvents));
650 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
651 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
652 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
653 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
655 // draw some example plots....
657 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
659 printf("CONTAINER NOT FOUND\n");
662 // projecting the containers to obtain histograms
663 // first argument = variable, second argument = step
666 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
667 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
668 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
669 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
670 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
671 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
672 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
673 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
674 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
675 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
676 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
677 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
679 // MC-Acceptance level
680 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
681 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
682 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
683 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
684 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
685 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
686 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
687 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
688 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
689 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
690 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
691 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
694 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
695 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
696 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
697 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
698 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
699 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
700 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
701 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
702 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
703 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
704 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
705 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
707 h00->SetTitle("pT_D0 (GeV/c)");
708 h10->SetTitle("rapidity");
709 h20->SetTitle("cosThetaStar");
710 h30->SetTitle("pT_pi (GeV/c)");
711 h40->SetTitle("pT_K (Gev/c)");
712 h50->SetTitle("cT (#mum)");
713 h60->SetTitle("dca (#mum)");
714 h70->SetTitle("d0_pi (#mum)");
715 h80->SetTitle("d0_K (#mum)");
716 h90->SetTitle("d0xd0 (#mum^2)");
717 h100->SetTitle("cosPointingAngle");
718 h100->SetTitle("phi (rad)");
720 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
721 h10->GetXaxis()->SetTitle("rapidity");
722 h20->GetXaxis()->SetTitle("cosThetaStar");
723 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
724 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
725 h50->GetXaxis()->SetTitle("cT (#mum)");
726 h60->GetXaxis()->SetTitle("dca (#mum)");
727 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
728 h80->GetXaxis()->SetTitle("d0_K (#mum)");
729 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
730 h100->GetXaxis()->SetTitle("cosPointingAngle");
731 h110->GetXaxis()->SetTitle("phi (rad)");
733 h01->SetTitle("pT_D0 (GeV/c)");
734 h11->SetTitle("rapidity");
735 h21->SetTitle("cosThetaStar");
736 h31->SetTitle("pT_pi (GeV/c)");
737 h41->SetTitle("pT_K (Gev/c)");
738 h51->SetTitle("cT (#mum)");
739 h61->SetTitle("dca (#mum)");
740 h71->SetTitle("d0_pi (#mum)");
741 h81->SetTitle("d0_K (#mum)");
742 h91->SetTitle("d0xd0 (#mum^2)");
743 h101->SetTitle("cosPointingAngle");
744 h111->GetXaxis()->SetTitle("phi (rad)");
746 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
747 h11->GetXaxis()->SetTitle("rapidity");
748 h21->GetXaxis()->SetTitle("cosThetaStar");
749 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
750 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
751 h51->GetXaxis()->SetTitle("cT (#mum)");
752 h61->GetXaxis()->SetTitle("dca (#mum)");
753 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
754 h81->GetXaxis()->SetTitle("d0_K (#mum)");
755 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
756 h101->GetXaxis()->SetTitle("cosPointingAngle");
757 h111->GetXaxis()->SetTitle("phi (rad)");
759 h04->SetTitle("pT_D0 (GeV/c)");
760 h14->SetTitle("rapidity");
761 h24->SetTitle("cosThetaStar");
762 h34->SetTitle("pT_pi (GeV/c)");
763 h44->SetTitle("pT_K (Gev/c)");
764 h54->SetTitle("cT (#mum)");
765 h64->SetTitle("dca (#mum)");
766 h74->SetTitle("d0_pi (#mum)");
767 h84->SetTitle("d0_K (#mum)");
768 h94->SetTitle("d0xd0 (#mum^2)");
769 h104->SetTitle("cosPointingAngle");
770 h114->GetXaxis()->SetTitle("phi (rad)");
772 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
773 h14->GetXaxis()->SetTitle("rapidity");
774 h24->GetXaxis()->SetTitle("cosThetaStar");
775 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
776 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
777 h54->GetXaxis()->SetTitle("cT (#mum)");
778 h64->GetXaxis()->SetTitle("dca (#mum)");
779 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
780 h84->GetXaxis()->SetTitle("d0_K (#mum)");
781 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
782 h104->GetXaxis()->SetTitle("cosPointingAngle");
783 h114->GetXaxis()->SetTitle("phi (rad)");
785 Double_t max0 = h00->GetMaximum();
786 Double_t max1 = h10->GetMaximum();
787 Double_t max2 = h20->GetMaximum();
788 Double_t max3 = h30->GetMaximum();
789 Double_t max4 = h40->GetMaximum();
790 Double_t max5 = h50->GetMaximum();
791 Double_t max6 = h60->GetMaximum();
792 Double_t max7 = h70->GetMaximum();
793 Double_t max8 = h80->GetMaximum();
794 Double_t max9 = h90->GetMaximum();
795 Double_t max10 = h100->GetMaximum();
796 Double_t max11 = h110->GetMaximum();
798 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
799 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
800 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
801 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
802 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
803 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
804 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
805 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
806 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
807 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
808 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
809 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
811 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
812 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
813 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
814 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
815 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
816 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
817 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
818 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
819 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
820 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
821 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
822 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
824 h00->SetMarkerStyle(20);
825 h10->SetMarkerStyle(24);
826 h20->SetMarkerStyle(21);
827 h30->SetMarkerStyle(25);
828 h40->SetMarkerStyle(27);
829 h50->SetMarkerStyle(28);
830 h60->SetMarkerStyle(20);
831 h70->SetMarkerStyle(24);
832 h80->SetMarkerStyle(21);
833 h90->SetMarkerStyle(25);
834 h100->SetMarkerStyle(27);
835 h110->SetMarkerStyle(28);
837 h00->SetMarkerColor(2);
838 h10->SetMarkerColor(2);
839 h20->SetMarkerColor(2);
840 h30->SetMarkerColor(2);
841 h40->SetMarkerColor(2);
842 h50->SetMarkerColor(2);
843 h60->SetMarkerColor(2);
844 h70->SetMarkerColor(2);
845 h80->SetMarkerColor(2);
846 h90->SetMarkerColor(2);
847 h100->SetMarkerColor(2);
848 h110->SetMarkerColor(2);
850 h01->SetMarkerStyle(20) ;
851 h11->SetMarkerStyle(24) ;
852 h21->SetMarkerStyle(21) ;
853 h31->SetMarkerStyle(25) ;
854 h41->SetMarkerStyle(27) ;
855 h51->SetMarkerStyle(28) ;
856 h61->SetMarkerStyle(20);
857 h71->SetMarkerStyle(24);
858 h81->SetMarkerStyle(21);
859 h91->SetMarkerStyle(25);
860 h101->SetMarkerStyle(27);
861 h111->SetMarkerStyle(28);
863 h01->SetMarkerColor(8);
864 h11->SetMarkerColor(8);
865 h21->SetMarkerColor(8);
866 h31->SetMarkerColor(8);
867 h41->SetMarkerColor(8);
868 h51->SetMarkerColor(8);
869 h61->SetMarkerColor(8);
870 h71->SetMarkerColor(8);
871 h81->SetMarkerColor(8);
872 h91->SetMarkerColor(8);
873 h101->SetMarkerColor(8);
874 h111->SetMarkerColor(8);
876 h04->SetMarkerStyle(20) ;
877 h14->SetMarkerStyle(24) ;
878 h24->SetMarkerStyle(21) ;
879 h34->SetMarkerStyle(25) ;
880 h44->SetMarkerStyle(27) ;
881 h54->SetMarkerStyle(28) ;
882 h64->SetMarkerStyle(20);
883 h74->SetMarkerStyle(24);
884 h84->SetMarkerStyle(21);
885 h94->SetMarkerStyle(25);
886 h104->SetMarkerStyle(27);
887 h114->SetMarkerStyle(28);
889 h04->SetMarkerColor(4);
890 h14->SetMarkerColor(4);
891 h24->SetMarkerColor(4);
892 h34->SetMarkerColor(4);
893 h44->SetMarkerColor(4);
894 h54->SetMarkerColor(4);
895 h64->SetMarkerColor(4);
896 h74->SetMarkerColor(4);
897 h84->SetMarkerColor(4);
898 h94->SetMarkerColor(4);
899 h104->SetMarkerColor(4);
900 h114->SetMarkerColor(4);
902 gStyle->SetCanvasColor(0);
903 gStyle->SetFrameFillColor(0);
904 gStyle->SetTitleFillColor(0);
905 gStyle->SetStatColor(0);
907 // drawing in 2 separate canvas for a matter of clearity
908 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
940 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
971 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1002 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1033 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1035 TH2D* corr1 =hcorr->Projection(0,2);
1036 TH2D* corr2 = hcorr->Projection(1,3);
1038 TCanvas * c7 =new TCanvas("c7","",800,400);
1041 corr1->Draw("text");
1043 corr2->Draw("text");
1046 TFile* file_projection = new TFile("file_projection.root","RECREATE");
1050 h00->Write("pT_D0_step0");
1051 h10->Write("rapidity_step0");
1052 h20->Write("cosThetaStar_step0");
1053 h30->Write("pT_pi_step0");
1054 h40->Write("pT_K_step0");
1055 h50->Write("cT_step0");
1056 h60->Write("dca_step0");
1057 h70->Write("d0_pi_step0");
1058 h80->Write("d0_K_step0");
1059 h90->Write("d0xd0_step0");
1060 h100->Write("cosPointingAngle_step0");
1061 h110->Write("phi_step0");
1063 h01->Write("pT_D0_step1");
1064 h11->Write("rapidity_step1");
1065 h21->Write("cosThetaStar_step1");
1066 h31->Write("pT_pi_step1");
1067 h41->Write("pT_K_step1");
1068 h51->Write("cT_step1");
1069 h61->Write("dca_step1");
1070 h71->Write("d0_pi_step1");
1071 h81->Write("d0_K_step1");
1072 h91->Write("d0xd0_step1");
1073 h101->Write("cosPointingAngle_step1");
1074 h111->Write("phi_step1");
1076 h04->Write("pT_D0_step2");
1077 h14->Write("rapidity_step2");
1078 h24->Write("cosThetaStar_step2");
1079 h34->Write("pT_pi_step2");
1080 h44->Write("pT_K_step2");
1081 h54->Write("cT_step2");
1082 h64->Write("dca_step2");
1083 h74->Write("d0_pi_step2");
1084 h80->Write("d0_K_step2");
1085 h94->Write("d0xd0_step2");
1086 h104->Write("cosPointingAngle_step2");
1087 h114->Write("phi_step2");
1089 file_projection->Close();
1094 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1095 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1096 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1097 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1099 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1100 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1101 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1102 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1106 //___________________________________________________________________________
1107 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1108 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1109 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1111 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1115 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
1118 //___________________________________________________________________________
1119 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1122 // to calculate cos(ThetaStar) for generated particle
1123 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1124 // (see where the function is called)
1127 Int_t pdgvtx = mcPart->GetPdgCode();
1128 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1129 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1130 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1131 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1132 AliDebug(2,"This is a D0");
1133 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1134 mcPartDaughter0 = mcPartDaughter1;
1135 mcPartDaughter1 = mcPartdummy;
1141 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1142 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1143 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1147 AliDebug(2,"D0bar");
1149 if (pdgprong0 == 211){
1150 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1151 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1152 mcPartDaughter0 = mcPartDaughter1;
1153 mcPartDaughter1 = mcPartdummy;
1154 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1155 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1158 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1159 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1161 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1162 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1164 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
1166 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1167 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1168 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1169 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1170 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1171 Double_t beta = p/e;
1172 Double_t gamma = e/massvtx;
1173 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1174 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1176 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1178 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1179 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1180 AliDebug(2,Form("cts = %f", cts));
1183 //___________________________________________________________________________
1184 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1187 // to calculate cT for generated particle
1190 Double_t xmcPart[3] = {0,0,0};
1191 Double_t xdaughter0[3] = {0,0,0};
1192 Double_t xdaughter1[3] = {0,0,0};
1193 mcPart->XvYvZv(xmcPart); // cm
1194 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1195 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1196 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1197 xmcPart[1]*xmcPart[1]+
1198 xmcPart[2]*xmcPart[2]);
1199 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1200 xdaughter0[1]*xdaughter0[1]+
1201 xdaughter0[2]*xdaughter0[2]);
1202 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1203 xdaughter1[1]*xdaughter1[1]+
1204 xdaughter1[2]*xdaughter1[2]);
1206 AliDebug(2, Form("D0: x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4));
1207 AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4));
1208 AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4));
1210 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1211 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1212 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1214 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1215 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1216 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1218 if ((cT0 - cT1)>1E-5) {
1219 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1222 // calculating cT from cT0
1224 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1225 Double_t p = mcPart-> P();
1226 Double_t cT = cT0*mass/p;
1227 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1228 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1229 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1232 //________________________________________________________________________________
1233 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1236 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1239 Bool_t isOk = kFALSE;
1241 // check whether the D0 decays in pi+K
1242 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1243 // could use a cut in the CF, but not clear how to define a TDecayChannel
1244 // implemented for the time being as a cut on the number of daughters - see later when
1245 // getting the daughters
1247 // getting the daughters
1248 Int_t daughter0 = mcPart->GetDaughter(0);
1249 Int_t daughter1 = mcPart->GetDaughter(1);
1250 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1251 if (daughter0 == 0 || daughter1 == 0) {
1252 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1255 if (TMath::Abs(daughter1 - daughter0) != 1) {
1256 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1259 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1260 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1261 if (!mcPartDaughter0 || !mcPartDaughter1) {
1262 AliWarning("At least one Daughter Particle not found in tree, skipping");
1265 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1266 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1267 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1268 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1269 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1273 Double_t vtx1[3] = {0,0,0}; // primary vertex
1274 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1275 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1276 mcPart->XvYvZv(vtx1); // cm
1277 // getting vertex from daughters
1278 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1279 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1280 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1281 AliError("Daughters have different secondary vertex, skipping the track");
1286 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1287 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1288 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1289 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1290 // inverting in case the positive daughter is the second one
1291 positiveDaugh = mcPartDaughter1;
1292 negativeDaugh = mcPartDaughter0;
1294 // getting the momentum from the daughters
1295 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1296 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1297 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1299 Double_t d0[2] = {0.,0.};
1301 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1303 Double_t cosThetaStar = 0.;
1304 Double_t cosThetaStarD0 = 0.;
1305 Double_t cosThetaStarD0bar = 0.;
1306 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1307 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1308 if (mcPart->GetPdgCode() == 421){ // D0
1309 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1310 cosThetaStar = cosThetaStarD0;
1312 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1313 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1314 cosThetaStar = cosThetaStarD0bar;
1317 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1320 if (cosThetaStar < -1 || cosThetaStar > 1) {
1321 AliWarning("Invalid value for cosine Theta star, returning");
1325 // calculate cos(Theta*) according to the method implemented herein
1327 Double_t cts = 9999.;
1328 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1329 if (cts < -1 || cts > 1) {
1330 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1332 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1333 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1336 Double_t cT = decay->Ct(421);
1338 // calculate cT from the method implemented herein
1340 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1342 if (TMath::Abs(cT1 - cT)>0.001) {
1343 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1346 // get the pT of the daughters
1351 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1352 pTpi = mcPartDaughter0->Pt();
1353 pTK = mcPartDaughter1->Pt();
1356 pTpi = mcPartDaughter1->Pt();
1357 pTK = mcPartDaughter0->Pt();
1360 vectorMC[0] = mcPart->Pt();
1361 vectorMC[1] = mcPart->Y() ;
1362 vectorMC[2] = cosThetaStar ;
1363 vectorMC[3] = pTpi ;
1365 vectorMC[5] = cT*1.E4 ; // in micron
1366 vectorMC[6] = mcPart->Phi() ;
1370 //_________________________________________________________________________________________________
1371 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1374 // checking whether the very mother of the D0 is a charm or a bottom quark
1377 Int_t pdgGranma = 0;
1379 mother = mcPart->GetMother();
1383 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1384 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1385 pdgGranma = mcGranma->GetPdgCode();
1386 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1387 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1388 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1391 mother = mcGranma->GetMother();