]> 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
Bug fix in PID step (ChiaraZ)
[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));
ae3d971c 652
653 // setting the use of the PID cut when applying the selection to FALSE - whatever it was. Keeping track of the original value
654 Bool_t iscutusingpid=fCuts->GetIsUsePID();
655 Int_t isselcuts=-1,isselpid=-1;
656 fCuts->SetUsePID(kFALSE);
ae3d971c 657 isselcuts = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kCandidate,aodEvent);
da576827 658 fCuts->SetUsePID(iscutusingpid); // restoring usage of the PID from the cuts object
ae3d971c 659 if (isselcuts == 3 || isselcuts == isD0D0bar){
65f0c9a8 660 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
ae3d971c 661 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR,fWeight) ;
f5ec6c30 662 icountRecoPPR++;
ae3e319c 663
664 if(!fAcceptanceUnf){ // unfolding
665
666 Double_t fill[4]; //fill response matrix
667
668 // dimensions 0&1 : pt,eta (Rec)
669
670 fill[0] = pt ;
671 fill[1] = rapidity;
672
673 // dimensions 2&3 : pt,eta (MC)
674
675 fill[2] = mcVtxHF->Pt();
676 fill[3] = mcVtxHF->Y();
677
678 fCorrelation->Fill(fill);
679
680 }
ae3d971c 681
ae3d971c 682 isselpid = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kPID);
683 if((fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==isD0D0bar)||(fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==3)){
174db2b5 684 AliDebug(2,"Particle passed PID cuts");
ae3d971c 685 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID,fWeight) ;
174db2b5 686 icountRecoPID++;
f5ec6c30 687 }
688 }
689 }
690 }
691 }
1e090d47 692 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
693 } // end loop on D0->Kpi
f5ec6c30 694
695 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
696
697 fCountReco+= icountReco;
ae3e319c 698 fCountVertex+= icountVertex;
699 fCountRefit+= icountRefit;
f5ec6c30 700 fCountRecoAcc+= icountRecoAcc;
701 fCountRecoITSClusters+= icountRecoITSClusters;
702 fCountRecoPPR+= icountRecoPPR;
174db2b5 703 fCountRecoPID+= icountRecoPID;
f5ec6c30 704
705 fHistEventsProcessed->Fill(0);
706 /* PostData(0) is taken care of by AliAnalysisTaskSE */
174db2b5 707 //PostData(1,fHistEventsProcessed) ;
708 //PostData(2,fCFManager->GetParticleContainer()) ;
709 //PostData(3,fCorrelation) ;
f5ec6c30 710}
711
712
713//___________________________________________________________________________
714void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
715{
716 // The Terminate() function is the last function to be called during
717 // a query. It always runs on the client, it can be used to present
718 // the results graphically or save the results to file.
719
f5ec6c30 720 AliAnalysisTaskSE::Terminate();
721
722 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
723 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
ae3e319c 724 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
725 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 726 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
727 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));
728 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));
729 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
730
731 // draw some example plots....
732
174db2b5 733 // AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
f5ec6c30 734 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
33aa234e 735 if(!cont) {
736 printf("CONTAINER NOT FOUND\n");
737 return;
738 }
f5ec6c30 739 // projecting the containers to obtain histograms
740 // first argument = variable, second argument = step
741
742 // MC-level
743 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
744 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
745 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
746 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
747 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
748 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
749 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
750 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
751 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
752 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
753 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
e8158e56 754 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
f5ec6c30 755
756 // MC-Acceptance level
757 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
758 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
759 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
760 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
761 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
762 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
763 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
764 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
765 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
766 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
767 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
e8158e56 768 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
f5ec6c30 769
770 // Reco-level
ae3e319c 771 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
772 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
773 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
774 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
775 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
776 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
777 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
778 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
779 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
780 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
781 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
782 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
f5ec6c30 783
784 h00->SetTitle("pT_D0 (GeV/c)");
785 h10->SetTitle("rapidity");
786 h20->SetTitle("cosThetaStar");
787 h30->SetTitle("pT_pi (GeV/c)");
788 h40->SetTitle("pT_K (Gev/c)");
789 h50->SetTitle("cT (#mum)");
790 h60->SetTitle("dca (#mum)");
791 h70->SetTitle("d0_pi (#mum)");
792 h80->SetTitle("d0_K (#mum)");
793 h90->SetTitle("d0xd0 (#mum^2)");
794 h100->SetTitle("cosPointingAngle");
e8158e56 795 h100->SetTitle("phi (rad)");
f5ec6c30 796
797 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
798 h10->GetXaxis()->SetTitle("rapidity");
799 h20->GetXaxis()->SetTitle("cosThetaStar");
800 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
801 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
802 h50->GetXaxis()->SetTitle("cT (#mum)");
803 h60->GetXaxis()->SetTitle("dca (#mum)");
804 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
805 h80->GetXaxis()->SetTitle("d0_K (#mum)");
806 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
807 h100->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 808 h110->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 809
810 h01->SetTitle("pT_D0 (GeV/c)");
811 h11->SetTitle("rapidity");
812 h21->SetTitle("cosThetaStar");
813 h31->SetTitle("pT_pi (GeV/c)");
814 h41->SetTitle("pT_K (Gev/c)");
815 h51->SetTitle("cT (#mum)");
816 h61->SetTitle("dca (#mum)");
817 h71->SetTitle("d0_pi (#mum)");
818 h81->SetTitle("d0_K (#mum)");
819 h91->SetTitle("d0xd0 (#mum^2)");
820 h101->SetTitle("cosPointingAngle");
e8158e56 821 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 822
823 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
824 h11->GetXaxis()->SetTitle("rapidity");
825 h21->GetXaxis()->SetTitle("cosThetaStar");
826 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
827 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
828 h51->GetXaxis()->SetTitle("cT (#mum)");
829 h61->GetXaxis()->SetTitle("dca (#mum)");
830 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
831 h81->GetXaxis()->SetTitle("d0_K (#mum)");
832 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
833 h101->GetXaxis()->SetTitle("cosPointingAngle");
e8158e56 834 h111->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 835
ae3e319c 836 h04->SetTitle("pT_D0 (GeV/c)");
837 h14->SetTitle("rapidity");
838 h24->SetTitle("cosThetaStar");
839 h34->SetTitle("pT_pi (GeV/c)");
840 h44->SetTitle("pT_K (Gev/c)");
841 h54->SetTitle("cT (#mum)");
842 h64->SetTitle("dca (#mum)");
843 h74->SetTitle("d0_pi (#mum)");
844 h84->SetTitle("d0_K (#mum)");
845 h94->SetTitle("d0xd0 (#mum^2)");
846 h104->SetTitle("cosPointingAngle");
847 h114->GetXaxis()->SetTitle("phi (rad)");
848
849 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
850 h14->GetXaxis()->SetTitle("rapidity");
851 h24->GetXaxis()->SetTitle("cosThetaStar");
852 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
853 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
854 h54->GetXaxis()->SetTitle("cT (#mum)");
855 h64->GetXaxis()->SetTitle("dca (#mum)");
856 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
857 h84->GetXaxis()->SetTitle("d0_K (#mum)");
858 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
859 h104->GetXaxis()->SetTitle("cosPointingAngle");
860 h114->GetXaxis()->SetTitle("phi (rad)");
f5ec6c30 861
862 Double_t max0 = h00->GetMaximum();
863 Double_t max1 = h10->GetMaximum();
864 Double_t max2 = h20->GetMaximum();
865 Double_t max3 = h30->GetMaximum();
866 Double_t max4 = h40->GetMaximum();
867 Double_t max5 = h50->GetMaximum();
868 Double_t max6 = h60->GetMaximum();
869 Double_t max7 = h70->GetMaximum();
870 Double_t max8 = h80->GetMaximum();
871 Double_t max9 = h90->GetMaximum();
872 Double_t max10 = h100->GetMaximum();
e8158e56 873 Double_t max11 = h110->GetMaximum();
f5ec6c30 874
875 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
876 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
877 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
878 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
879 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
880 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
881 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
882 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
883 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
884 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
885 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 886 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 887
888 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
889 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
890 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
891 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
892 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
893 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
894 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
895 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
896 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
897 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
898 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
e8158e56 899 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
f5ec6c30 900
901 h00->SetMarkerStyle(20);
902 h10->SetMarkerStyle(24);
903 h20->SetMarkerStyle(21);
904 h30->SetMarkerStyle(25);
905 h40->SetMarkerStyle(27);
906 h50->SetMarkerStyle(28);
907 h60->SetMarkerStyle(20);
908 h70->SetMarkerStyle(24);
909 h80->SetMarkerStyle(21);
910 h90->SetMarkerStyle(25);
911 h100->SetMarkerStyle(27);
e8158e56 912 h110->SetMarkerStyle(28);
f5ec6c30 913
914 h00->SetMarkerColor(2);
915 h10->SetMarkerColor(2);
916 h20->SetMarkerColor(2);
917 h30->SetMarkerColor(2);
918 h40->SetMarkerColor(2);
919 h50->SetMarkerColor(2);
920 h60->SetMarkerColor(2);
921 h70->SetMarkerColor(2);
922 h80->SetMarkerColor(2);
923 h90->SetMarkerColor(2);
924 h100->SetMarkerColor(2);
e8158e56 925 h110->SetMarkerColor(2);
f5ec6c30 926
927 h01->SetMarkerStyle(20) ;
928 h11->SetMarkerStyle(24) ;
929 h21->SetMarkerStyle(21) ;
930 h31->SetMarkerStyle(25) ;
931 h41->SetMarkerStyle(27) ;
932 h51->SetMarkerStyle(28) ;
933 h61->SetMarkerStyle(20);
934 h71->SetMarkerStyle(24);
935 h81->SetMarkerStyle(21);
936 h91->SetMarkerStyle(25);
937 h101->SetMarkerStyle(27);
e8158e56 938 h111->SetMarkerStyle(28);
f5ec6c30 939
940 h01->SetMarkerColor(8);
941 h11->SetMarkerColor(8);
942 h21->SetMarkerColor(8);
943 h31->SetMarkerColor(8);
944 h41->SetMarkerColor(8);
945 h51->SetMarkerColor(8);
946 h61->SetMarkerColor(8);
947 h71->SetMarkerColor(8);
948 h81->SetMarkerColor(8);
949 h91->SetMarkerColor(8);
950 h101->SetMarkerColor(8);
e8158e56 951 h111->SetMarkerColor(8);
f5ec6c30 952
ae3e319c 953 h04->SetMarkerStyle(20) ;
954 h14->SetMarkerStyle(24) ;
955 h24->SetMarkerStyle(21) ;
956 h34->SetMarkerStyle(25) ;
957 h44->SetMarkerStyle(27) ;
958 h54->SetMarkerStyle(28) ;
959 h64->SetMarkerStyle(20);
960 h74->SetMarkerStyle(24);
961 h84->SetMarkerStyle(21);
962 h94->SetMarkerStyle(25);
963 h104->SetMarkerStyle(27);
964 h114->SetMarkerStyle(28);
965
966 h04->SetMarkerColor(4);
967 h14->SetMarkerColor(4);
968 h24->SetMarkerColor(4);
969 h34->SetMarkerColor(4);
970 h44->SetMarkerColor(4);
971 h54->SetMarkerColor(4);
972 h64->SetMarkerColor(4);
973 h74->SetMarkerColor(4);
974 h84->SetMarkerColor(4);
975 h94->SetMarkerColor(4);
976 h104->SetMarkerColor(4);
977 h114->SetMarkerColor(4);
f5ec6c30 978
979 gStyle->SetCanvasColor(0);
980 gStyle->SetFrameFillColor(0);
981 gStyle->SetTitleFillColor(0);
982 gStyle->SetStatColor(0);
983
984 // drawing in 2 separate canvas for a matter of clearity
985 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
986 c1->Divide(3,3);
987
988 c1->cd(1);
989 h00->Draw("p");
990 c1->cd(1);
991 c1->cd(2);
992 h01->Draw("p");
993 c1->cd(2);
994 c1->cd(3);
ae3e319c 995 h04->Draw("p");
f5ec6c30 996 c1->cd(3);
997 c1->cd(4);
998 h10->Draw("p");
999 c1->cd(4);
1000 c1->cd(5);
1001 h11->Draw("p");
1002 c1->cd(5);
1003 c1->cd(6);
ae3e319c 1004 h14->Draw("p");
f5ec6c30 1005 c1->cd(6);
1006 c1->cd(7);
1007 h20->Draw("p");
1008 c1->cd(7);
1009 c1->cd(8);
1010 h21->Draw("p");
1011 c1->cd(8);
1012 c1->cd(9);
ae3e319c 1013 h24->Draw("p");
f5ec6c30 1014 c1->cd(9);
1015 c1->cd();
1016
1017 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1018 c2->Divide(3,3);
1019 c2->cd(1);
1020 h30->Draw("p");
1021 c2->cd(1);
1022 c2->cd(2);
1023 h31->Draw("p");
1024 c2->cd(2);
1025 c2->cd(3);
ae3e319c 1026 h34->Draw("p");
f5ec6c30 1027 c2->cd(3);
1028 c2->cd(4);
1029 h40->Draw("p");
1030 c2->cd(4);
1031 c2->cd(5);
1032 h41->Draw("p");
1033 c2->cd(5);
1034 c2->cd(6);
ae3e319c 1035 h44->Draw("p");
f5ec6c30 1036 c2->cd(6);
1037 c2->cd(7);
1038 h50->Draw("p");
1039 c2->cd(7);
1040 c2->cd(8);
1041 h51->Draw("p");
1042 c2->cd(8);
1043 c2->cd(9);
ae3e319c 1044 h54->Draw("p");
f5ec6c30 1045 c2->cd(9);
1046 c2->cd();
1047
1048 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1049 c3->Divide(3,3);
1050 c3->cd(1);
1051 h60->Draw("p");
1052 c3->cd(1);
1053 c3->cd(2);
1054 h61->Draw("p");
1055 c3->cd(2);
1056 c3->cd(3);
ae3e319c 1057 h64->Draw("p");
f5ec6c30 1058 c3->cd(3);
1059 c3->cd(4);
1060 h70->Draw("p");
1061 c3->cd(4);
1062 c3->cd(5);
1063 h71->Draw("p");
1064 c3->cd(5);
1065 c3->cd(6);
ae3e319c 1066 h74->Draw("p");
f5ec6c30 1067 c3->cd(6);
1068 c3->cd(7);
1069 h80->Draw("p");
1070 c3->cd(7);
1071 c3->cd(8);
1072 h81->Draw("p");
1073 c3->cd(8);
1074 c3->cd(9);
ae3e319c 1075 h84->Draw("p");
f5ec6c30 1076 c3->cd(9);
1077 c3->cd();
1078
e8158e56 1079 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1080 c4->Divide(3,3);
f5ec6c30 1081 c4->cd(1);
1082 h90->Draw("p");
1083 c4->cd(1);
1084 c4->cd(2);
1085 h91->Draw("p");
1086 c4->cd(2);
1087 c4->cd(3);
ae3e319c 1088 h94->Draw("p");
f5ec6c30 1089 c4->cd(3);
1090 c4->cd(4);
1091 h100->Draw("p");
1092 c4->cd(4);
1093 c4->cd(5);
1094 h101->Draw("p");
1095 c4->cd(5);
1096 c4->cd(6);
ae3e319c 1097 h104->Draw("p");
f5ec6c30 1098 c4->cd(6);
e8158e56 1099 c4->cd(7);
1100 h110->Draw("p");
1101 c4->cd(7);
1102 c4->cd(8);
1103 h111->Draw("p");
1104 c4->cd(8);
1105 c4->cd(9);
ae3e319c 1106 h114->Draw("p");
e8158e56 1107 c4->cd(9);
f5ec6c30 1108 c4->cd();
1109
a5d92cb4 1110 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1111
1112 TH2D* corr1 =hcorr->Projection(0,2);
1113 TH2D* corr2 = hcorr->Projection(1,3);
1114
1115 TCanvas * c7 =new TCanvas("c7","",800,400);
1116 c7->Divide(2,1);
1117 c7->cd(1);
1118 corr1->Draw("text");
1119 c7->cd(2);
1120 corr2->Draw("text");
1121
1122
c59e1898 1123 TString projection_filename="CFtaskHFprojection";
1124 if(fKeepD0fromBOnly) projection_filename+="_KeepD0fromBOnly";
1125 else if(fKeepD0fromB) projection_filename+="_KeepD0fromB";
1126 projection_filename+=".root";
1127 TFile* file_projection = new TFile(projection_filename.Data(),"RECREATE");
a5d92cb4 1128
1129 corr1->Write();
1130 corr2->Write();
844d235c 1131 h00->Write("pT_D0_step0");
1132 h10->Write("rapidity_step0");
1133 h20->Write("cosThetaStar_step0");
1134 h30->Write("pT_pi_step0");
1135 h40->Write("pT_K_step0");
1136 h50->Write("cT_step0");
1137 h60->Write("dca_step0");
1138 h70->Write("d0_pi_step0");
1139 h80->Write("d0_K_step0");
1140 h90->Write("d0xd0_step0");
1141 h100->Write("cosPointingAngle_step0");
e8158e56 1142 h110->Write("phi_step0");
844d235c 1143
1144 h01->Write("pT_D0_step1");
1145 h11->Write("rapidity_step1");
1146 h21->Write("cosThetaStar_step1");
1147 h31->Write("pT_pi_step1");
1148 h41->Write("pT_K_step1");
1149 h51->Write("cT_step1");
1150 h61->Write("dca_step1");
1151 h71->Write("d0_pi_step1");
1152 h81->Write("d0_K_step1");
1153 h91->Write("d0xd0_step1");
1154 h101->Write("cosPointingAngle_step1");
e8158e56 1155 h111->Write("phi_step1");
844d235c 1156
ae3e319c 1157 h04->Write("pT_D0_step2");
1158 h14->Write("rapidity_step2");
1159 h24->Write("cosThetaStar_step2");
1160 h34->Write("pT_pi_step2");
1161 h44->Write("pT_K_step2");
1162 h54->Write("cT_step2");
1163 h64->Write("dca_step2");
1164 h74->Write("d0_pi_step2");
844d235c 1165 h80->Write("d0_K_step2");
ae3e319c 1166 h94->Write("d0xd0_step2");
1167 h104->Write("cosPointingAngle_step2");
1168 h114->Write("phi_step2");
844d235c 1169
1170 file_projection->Close();
1171
a5d92cb4 1172
1173
f5ec6c30 1174 /*
1175 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1176 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1177 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1178 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1179
1180 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1181 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1182 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1183 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
a8b9a95a 1184 */
f5ec6c30 1185}
1186
1187//___________________________________________________________________________
1188void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1189 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1190 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1191 //
1192 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1193
1194 //slot #1
1195 OpenFile(1);
a8b9a95a 1196 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
f5ec6c30 1197}
1198
1199//___________________________________________________________________________
1200Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1201
1202 //
1203 // to calculate cos(ThetaStar) for generated particle
1204 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1205 // (see where the function is called)
1206 //
1207
1208 Int_t pdgvtx = mcPart->GetPdgCode();
1209 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1210 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1211 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1212 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1213 AliDebug(2,"This is a D0");
1214 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1215 mcPartDaughter0 = mcPartDaughter1;
1216 mcPartDaughter1 = mcPartdummy;
1217 }
1218 else{
1219 AliInfo("D0bar");
1220 }
1221 */
1222 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1223 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1224 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1225 AliDebug(2,"D0");
1226 }
1227 else{
1228 AliDebug(2,"D0bar");
1229 }
1230 if (pdgprong0 == 211){
1231 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1232 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1233 mcPartDaughter0 = mcPartDaughter1;
1234 mcPartDaughter1 = mcPartdummy;
1235 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1236 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1237 }
1238
1239 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1240 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1241 Double_t massp[2];
1242 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1243 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1244
1245 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);
1246
1247 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1248 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1249 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1250 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1251 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1252 Double_t beta = p/e;
1253 Double_t gamma = e/massvtx;
1254 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1255 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1256
1257 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1258
1259 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1260 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1261 AliDebug(2,Form("cts = %f", cts));
1262 return cts;
1263}
1264//___________________________________________________________________________
1265Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1266
1267 //
1268 // to calculate cT for generated particle
1269 //
1270
1271 Double_t xmcPart[3] = {0,0,0};
1272 Double_t xdaughter0[3] = {0,0,0};
1273 Double_t xdaughter1[3] = {0,0,0};
1274 mcPart->XvYvZv(xmcPart); // cm
1275 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1276 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1277 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1278 xmcPart[1]*xmcPart[1]+
1279 xmcPart[2]*xmcPart[2]);
1280 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1281 xdaughter0[1]*xdaughter0[1]+
1282 xdaughter0[2]*xdaughter0[2]);
1283 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1284 xdaughter1[1]*xdaughter1[1]+
1285 xdaughter1[2]*xdaughter1[2]);
1286
1287 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));
1288 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));
1289 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));
1290
1291 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1292 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1293 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1294
1295 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1296 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1297 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1298
e8158e56 1299 if ((cT0 - cT1)>1E-5) {
1300 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 1301 }
1302
1303 // calculating cT from cT0
1304
1305 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1306 Double_t p = mcPart-> P();
1307 Double_t cT = cT0*mass/p;
1308 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1309 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1310 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1311 return cT;
1312}
1313//________________________________________________________________________________
1314Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1315
1316 //
1317 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1318 //
1319
1320 Bool_t isOk = kFALSE;
1321
1322 // check whether the D0 decays in pi+K
1323 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1324 // could use a cut in the CF, but not clear how to define a TDecayChannel
1325 // implemented for the time being as a cut on the number of daughters - see later when
1326 // getting the daughters
1327
1328 // getting the daughters
1329 Int_t daughter0 = mcPart->GetDaughter(0);
1330 Int_t daughter1 = mcPart->GetDaughter(1);
1331 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1332 if (daughter0 == 0 || daughter1 == 0) {
1333 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1334 return isOk;
1335 }
421623bd 1336 if (TMath::Abs(daughter1 - daughter0) != 1) {
f5ec6c30 1337 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1338 return isOk;
1339 }
1340 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1341 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1342 if (!mcPartDaughter0 || !mcPartDaughter1) {
1343 AliWarning("At least one Daughter Particle not found in tree, skipping");
1344 return isOk;
1345 }
2bcc508a 1346 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1347 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1348 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1349 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1350 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1351 return isOk;
1352 }
1353
f5ec6c30 1354 Double_t vtx1[3] = {0,0,0}; // primary vertex
1355 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1356 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1357 mcPart->XvYvZv(vtx1); // cm
1358 // getting vertex from daughters
1359 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1360 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1361 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1362 AliError("Daughters have different secondary vertex, skipping the track");
1363 return isOk;
1364 }
1365 Int_t nprongs = 2;
1366 Short_t charge = 0;
1367 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1368 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1369 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1370 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1371 // inverting in case the positive daughter is the second one
1372 positiveDaugh = mcPartDaughter1;
1373 negativeDaugh = mcPartDaughter0;
1374 }
1375 // getting the momentum from the daughters
1376 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1377 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1378 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1379
1380 Double_t d0[2] = {0.,0.};
1381
1382 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1383
1384 Double_t cosThetaStar = 0.;
1385 Double_t cosThetaStarD0 = 0.;
1386 Double_t cosThetaStarD0bar = 0.;
1387 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1388 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1389 if (mcPart->GetPdgCode() == 421){ // D0
1390 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1391 cosThetaStar = cosThetaStarD0;
1392 }
1393 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1394 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1395 cosThetaStar = cosThetaStarD0bar;
1396 }
1397 else{
1398 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1399 return vectorMC;
1400 }
1401 if (cosThetaStar < -1 || cosThetaStar > 1) {
1402 AliWarning("Invalid value for cosine Theta star, returning");
1403 return isOk;
1404 }
1405
1406 // calculate cos(Theta*) according to the method implemented herein
1407
1408 Double_t cts = 9999.;
1409 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1410 if (cts < -1 || cts > 1) {
1411 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1412 }
1413 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1414 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1415 }
1416
1417 Double_t cT = decay->Ct(421);
1418
1419 // calculate cT from the method implemented herein
1420 Double_t cT1 = 0.;
1421 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1422
1423 if (TMath::Abs(cT1 - cT)>0.001) {
1424 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1425 }
1426
1427 // get the pT of the daughters
1428
1429 Double_t pTpi = 0.;
1430 Double_t pTK = 0.;
1431
1432 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1433 pTpi = mcPartDaughter0->Pt();
1434 pTK = mcPartDaughter1->Pt();
1435 }
1436 else {
1437 pTpi = mcPartDaughter1->Pt();
1438 pTK = mcPartDaughter0->Pt();
1439 }
1440
1441 vectorMC[0] = mcPart->Pt();
1442 vectorMC[1] = mcPart->Y() ;
1443 vectorMC[2] = cosThetaStar ;
1444 vectorMC[3] = pTpi ;
1445 vectorMC[4] = pTK ;
1446 vectorMC[5] = cT*1.E4 ; // in micron
e8158e56 1447 vectorMC[6] = mcPart->Phi() ;
f5ec6c30 1448 isOk = kTRUE;
1449 return isOk;
1450}
e8158e56 1451//_________________________________________________________________________________________________
1452Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1453
1454 //
1455 // checking whether the very mother of the D0 is a charm or a bottom quark
1456 //
1457
1458 Int_t pdgGranma = 0;
1459 Int_t mother = 0;
1460 mother = mcPart->GetMother();
1461 Int_t istep = 0;
1462 while (mother >0 ){
1463 istep++;
1464 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1465 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1466 pdgGranma = mcGranma->GetPdgCode();
1467 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1468 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1469 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1470 break;
1471 }
1472 mother = mcGranma->GetMother();
1473 }
1474 return pdgGranma;
1475}
ae3d971c 1476//__________________________________________________________________________________________________
1477Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt){
1478
1479 //
1480 // calculating the weight to fill the container
1481 //
1482
1483 // FNOLL central:
1484 // p0 = 0.2149
1485 // p1 = 2.433
1486 // p2 = 2.102
1487 // p3 = 2.5
1488
1489 // FNOLL min:
1490 // p0 = 0.06947
1491 // p1 = 4.795
1492 // p2 = 2.878
1493 // p3 = 2.5
1494
1495 // FNOLL max:
1496 // p0 = 0.2739
1497 // p1 = 2.136
1498 // p2 = 2.066
1499 // p3 = 2.5
1500
1501 // PYTHIA
1502 // p0 = 3.63832E-2
1503 // p1 = 1.88238
1504 // p2 = 1.34854
1505 // p3 = 2.5
1506
1507 Double_t func1[4] = {0.2149,2.433,2.102,2.5};
1508 Double_t func2[4] = {3.63832E-2,1.88238,1.34854,2.5};
1509
1510 Double_t dndpt_func1 = dNdptFit(pt,func1);
1511 Double_t dndpt_func2 = dNdptFit(pt,func2);
1512 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1513 return dndpt_func1/dndpt_func2;
1514}
1515
1516//__________________________________________________________________________________________________
1517Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::dNdptFit(Float_t pt, Double_t* par){
1518
1519 //
1520 // calculating dNdpt
1521 //
1522
1523 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1524 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1525
1526 return dNdpt;
1527}