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