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