AliPhysics  64f4410 (64f4410)
AliCFTaskForDStarAnalysis.cxx
Go to the documentation of this file.
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 /* $Id$ */
17 
18 //-----------------------------------------------------------------------
19 // Class for DStar corrections:
20 //
21 // The D0 cutting varibles position in the container and container
22 // binning method from a C.Zampolli example
23 // In this way a simple comparison between D0 and D0 from D* is possible.
24 //
25 //-----------------------------------------------------------------------
26 // Author : A.Grelli, Utrecht University
27 //
28 // a.grelli@uu.nl
29 //-----------------------------------------------------------------------
30 
31 #include <TCanvas.h>
32 #include <TParticle.h>
33 #include <TDatabasePDG.h>
34 #include <TH1I.h>
35 #include <TStyle.h>
36 #include <TFile.h>
37 
39 #include "AliMCEvent.h"
40 #include "AliCFManager.h"
41 #include "AliCFContainer.h"
42 #include "AliLog.h"
43 #include "AliAnalysisManager.h"
44 #include "AliAODHandler.h"
45 #include "AliAODEvent.h"
46 #include "AliAODRecoDecay.h"
47 #include "AliAODRecoDecayHF.h"
48 #include "AliAODRecoCascadeHF.h"
50 #include "AliAODMCParticle.h"
51 #include "AliAODMCHeader.h"
52 #include "AliESDtrack.h"
53 #include "TChain.h"
54 #include "THnSparse.h"
55 #include "TH2D.h"
56 
57 //__________________________________________________________________________
60  fCFManager(0x0),
61  fHistEventsProcessed(0x0),
62  fCorrelation(0x0),
63  fCountRecoDStarSel(0),
64  fEvents(0),
65  fMinITSClusters(5),
66  fMinITSClustersSoft(4),
67  fAcceptanceUnf(kTRUE)
68 {
69  //
70  //Default ctor
71  //
72 }
73 //___________________________________________________________________________
75  AliAnalysisTaskSE(name),
76  fCFManager(0x0),
78  fCorrelation(0x0),
80  fEvents(0),
81  fMinITSClusters(5),
83  fAcceptanceUnf(kTRUE)
84 {
85  //
86  // Constructor. Initialization of Inputs and Outputs
87  //
88  Info("AliCFTaskForDStarAnalysis","Calling Constructor");
89 
90  DefineOutput(1,TH1I::Class());
91  DefineOutput(2,AliCFContainer::Class());
92  DefineOutput(3,THnSparseD::Class());
93 }
94 
95 //___________________________________________________________________________
97 {
98  //
99  // Assignment operator
100  //
101  if (this!=&c) {
102  AliAnalysisTaskSE::operator=(c) ;
105  }
106  return *this;
107 }
108 
109 //___________________________________________________________________________
116  fEvents(c.fEvents),
120 {
121  //
122  // Copy Constructor
123  //
124 }
125 
126 //___________________________________________________________________________
128  //
129  //destructor
130  //
131  if (fCFManager) delete fCFManager ;
133  if (fCorrelation) delete fCorrelation ;
134 }
135 
136 //_________________________________________________
138 {
139  //
140  // Main loop function
141  //
142 
143  if (!fInputEvent) {
144  Error("UserExec","NO EVENT FOUND!");
145  return;
146  }
147 
148  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
149 
150  TClonesArray *arrayDStartoD0pi=0;
151 
152  if(!aodEvent && AODEvent() && IsStandardAOD()) {
153  // In case there is an AOD handler writing a standard AOD, use the AOD
154  // event in memory rather than the input (ESD) event.
155  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
156  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
157  // have to taken from the AOD event hold by the AliAODExtension
158  AliAODHandler* aodHandler = (AliAODHandler*)
159  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
160  if(aodHandler->GetExtensions()) {
161  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
162  AliAODEvent *aodFromExt = ext->GetAOD();
163  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
164  }
165  } else {
166  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
167  }
168 
169  if (!arrayDStartoD0pi) {
170  AliError("Could not find array of HF vertices");
171  return;
172  }
173 
174  // fix for temporary bug in ESDfilter
175  // the AODs with null vertex pointer didn't pass the PhysSel
176  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
177 
178  fEvents++;
179  if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
180 
181  fCFManager->SetRecEventInfo(aodEvent);
182  fCFManager->SetMCEventInfo(aodEvent);
183 
184  // event selection
185  Double_t containerInput[15] ;
186  Double_t containerInputMC[15] ;
187 
188  //loop on the MC event
189 
190  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
191  if (!mcArray) {
192  AliError("Could not find Monte-Carlo in AOD");
193  return;
194  }
195 
196  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
197  if (!mcHeader) {
198  AliError("Could not find MC Header in AOD");
199  return;
200  }
201 
202  // AOD primary vertex
203  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
204  Double_t zPrimVertex = vtx1->GetZ();
205  Double_t zMCVertex = mcHeader->GetVtxZ();
206  Bool_t vtxFlag = kTRUE;
207  TString title=vtx1->GetTitle();
208  if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
209 
210  for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
211  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
212  if (!mcPart) {
213  AliWarning("MC Particle not found in tree, skipping");
214  continue;
215  }
216 
217  // check the MC-level cuts
218  if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
219 
220  // fill the container for Gen-level selection
221  Double_t vectorMC[9] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
222 
223  if (GetDStarMCParticle(mcPart, mcArray, vectorMC)){
224 
225  containerInputMC[0] = vectorMC[0];
226  containerInputMC[1] = vectorMC[1] ;
227  containerInputMC[2] = vectorMC[2] ;
228  containerInputMC[3] = vectorMC[3] ;
229  containerInputMC[4] = vectorMC[4] ;
230  containerInputMC[5] = vectorMC[5] ;
231  containerInputMC[6] = 0.;
232  containerInputMC[7] = 0.;
233  containerInputMC[8] = 0.;
234  containerInputMC[9] = -100000.;
235  containerInputMC[10] = 1.01;
236  containerInputMC[11] = vectorMC[6];
237  containerInputMC[12] = zMCVertex; // z vertex
238  containerInputMC[13] = vectorMC[7]; // pt D0 pion
239  containerInputMC[14] = vectorMC[8]; // pt D0 kaon
240 
241  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
242 
243  // check the MC-Acceptance level cuts
244  // since standard CF functions are not applicable, using Kine Cuts on daughters
245 
246  Int_t daughter0 = mcPart->GetDaughter(0);
247  Int_t daughter1 = mcPart->GetDaughter(1);
248 
249  AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
250 
251  //D0
252  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
253  if (!mcPartDaughter0) {
254  AliError("Could not find Monte-Carlo in AOD");
255  return;
256  }
257 
258  // Soft Pion
259  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
260  if (!mcPartDaughter1) {
261  AliError("Could not find Monte-Carlo in AOD");
262  return;
263  }
264 
265  // Acceptance variables for the soft pion
266  Double_t eta1 = mcPartDaughter1->Eta();
267  Double_t pt1 = mcPartDaughter1->Pt();
268 
269  Int_t daughD00 = 0;
270  Int_t daughD01 = 0;
271 
272  // Just to be sure to take the right particles
273  if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
274  daughD00 = mcPartDaughter0->GetDaughter(0);
275  daughD01 = mcPartDaughter0->GetDaughter(1);
276  }else{
277  daughD00 = mcPartDaughter1->GetDaughter(0);
278  daughD01 = mcPartDaughter1->GetDaughter(1);
279  }
280 
281  // the D0 daughters
282  AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
283  AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
284 
285  if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
286  AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
287  return;
288  }
289 
290  // D0 daughters - needed for acceptance
291  Double_t theD0pt0 = mcD0PartDaughter0->Pt();
292  Double_t theD0pt1 = mcD0PartDaughter1->Pt();
293  Double_t theD0eta0 = mcD0PartDaughter0->Eta();
294  Double_t theD0eta1 = mcD0PartDaughter1->Eta();
295 
296  // ACCEPTANCE REQUESTS ---------
297 
298  // soft pion
299  Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
300  // Do daughters
301  Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1);
302  Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
303 
304  if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
305 
306  AliDebug(2, "D* Daughter particles in acceptance");
307 
308  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
309 
310  // check on the vertex
311  if (vtxFlag){
312  // filling the container if the vertex is ok
313  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
314 
315  Bool_t refitFlag = kTRUE;
316  for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
317  AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
318 
319  // refit only for D0 daughters
320  if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
321  if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
322  refitFlag = kFALSE;
323  }
324  }
325  }
326  if (refitFlag){ // refit
327 
328  fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
329 
330  } // end of refit
331  } // end of vertex
332  } //end of acceptance
333  } // end of MC D*
334  else {
335  AliDebug(3,"Problems in filling the container");
336  continue;
337  }
338  } // end of MC loop
339 
340  //rec
341  AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
342 
343  //D* and D0 prongs needed to MatchToMC method
344  Int_t pdgDgDStartoD0pi[2]={421,211};
345  Int_t pdgDgD0toKpi[2]={321,211};
346 
347  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
348 
349  // D* candidates
350  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
351 
352  // D0 from the reco cascade
353  AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
354  Bool_t unsetvtx=kFALSE;
355 
356  // needed for pointing angle
357  if(!D0particle->GetOwnPrimaryVtx()) {
358  D0particle->SetOwnPrimaryVtx(vtx1);
359  unsetvtx=kTRUE;
360  }
361 
362  // find associated MC particle for D* ->D0toKpi
363  Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
364 
365  // find D0->Kpi ... needed in the following
366  Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
367 
368  if (mcLabel == -1 || mcD0Label ==-1)
369  {
370  AliDebug(2,"No MC particle found");
371  }
372  else {
373 
374  // the D* and the D0 in MC
375  AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
376  AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
377 
378  if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
379  AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
380  continue;
381  }
382 
383  // soft pion
384  AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
385 
386  //D0tokpi
387  AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
388  AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
389 
390  // check if associated MC v0 passes the cuts
391  if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) {
392  AliDebug(2, "Skipping the particles due to cuts");
393  continue;
394  }
395 
396  // fill the container...
397  Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
398  Double_t rapidity = dstarD0pi->YDstar();
399  Double_t cosThetaStar = 9999.;
400  Double_t pTpi = 0.;
401  Double_t pTK = 0.;
402  Double_t dca = (D0particle->GetDCA())*1E4;
403  Double_t d0pi = 0.;
404  Double_t d0K = 0.;
405  Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
406  Double_t cosPointingAngle = D0particle->CosPointingAngle();
407  Double_t phi = D0particle->Phi();
408 
409  // Select D0 cutting variables
410  Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
411 
412  // D0 related variables
413  if (pdgCode > 0){
414 
415  cosThetaStar = D0particle->CosThetaStarD0();
416  pTpi = D0particle->PtProng(0);
417  pTK = D0particle->PtProng(1);
418  d0pi = (D0particle->Getd0Prong(0))*1E4;
419  d0K = (D0particle->Getd0Prong(1))*1E4;
420 
421  }
422  else {
423 
424  cosThetaStar = D0particle->CosThetaStarD0bar();
425  pTpi = D0particle->PtProng(1);
426  pTK = D0particle->PtProng(0);
427  d0pi = (D0particle->Getd0Prong(1))*1E4;
428  d0K = (D0particle->Getd0Prong(0))*1E4;
429 
430  }
431 
432  // ct of the D0 from D*
433  Double_t cT = D0particle->CtD0();
434 
435  containerInput[0] = pt;
436  containerInput[1] = rapidity;
437  containerInput[2] = cosThetaStar;
438  containerInput[3] = track2->Pt();
439  containerInput[4] = D0particle->Pt();
440  containerInput[5] = cT*1.E4; // in micron
441  containerInput[6] = dca; // in micron
442  containerInput[7] = d0pi; // in micron
443  containerInput[8] = d0K; // in micron
444  containerInput[9] = d0xd0; // in micron^2
445  containerInput[10] = cosPointingAngle; // in micron
446  containerInput[11] = phi;
447  containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
448  containerInput[13] = pTpi; // D0 pion
449  containerInput[14] = pTK; // D0 kaon
450 
451  fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
452 
453  // refit in ITS and TPC for D0 daughters
454  if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
455  continue;
456  }
457 
458  // reft in ITS for soft pion
459  if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
460  continue;
461  }
462 
463  // cut in acceptance for the soft pion and for the D0 daughters
464  Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
465  Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
466 
467  // soft pion acceptance ... is it fine 0.9?????
468  Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
469 
470  if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
471  AliDebug(2,"D* reco daughters in acceptance");
472  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
473 
474  if(fAcceptanceUnf){
475  Double_t fill[4]; //fill response matrix
476 
477  // dimensions 0&1 : pt,eta (Rec)
478  fill[0] = pt ;
479  fill[1] = rapidity;
480  // dimensions 2&3 : pt,eta (MC)
481  fill[2] = mcVtxHFDStar->Pt();
482  fill[3] = mcVtxHFDStar->Y();
483  fCorrelation->Fill(fill);
484  }
485 
486  // cut on the min n. of clusters in ITS for the D0 and soft pion
487  Int_t ncls0=0,ncls1=0,ncls2=0;
488  for(Int_t l=0;l<6;l++) {
489  if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
490  if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
491  if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
492  }
493  // see AddTask for soft pion ITS clusters request
494  AliDebug(2, Form("n clusters = %d", ncls0));
495 
496  if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
497  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
498 
499  // D0 cuts optimized for D* analysis
500  Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
501 
502  // needed for cuts
503  Double_t theD0pt = D0particle->Pt();
504 
505  if (theD0pt <= 1){ // first bin not optimized
506  cuts[0] = 400;
507  cuts[1] = 0.8;
508  cuts[2] = 0.21;
509  cuts[3] = 500;
510  cuts[4] = 500;
511  cuts[5] = -20000;
512  cuts[6] = 0.6;
513  }
514  else if (theD0pt > 1 && theD0pt <= 2){
515  cuts[0] = 200;
516  cuts[1] = 0.7;
517  cuts[2] = 0.8;
518  cuts[3] = 210;
519  cuts[4] = 210;
520  cuts[5] = -20000;
521  cuts[6] = 0.9;
522  }
523  else if (theD0pt > 2 && theD0pt <= 3){
524  cuts[0] = 400;
525  cuts[1] = 0.8;
526  cuts[2] = 0.8;
527  cuts[3] = 420;
528  cuts[4] = 350;
529  cuts[5] = -8500;
530  cuts[6] = 0.9;
531  }
532  else if (theD0pt > 3 && theD0pt <= 5){
533  cuts[0] = 160;
534  cuts[1] = 1.0;
535  cuts[2] = 1.2;
536  cuts[3] = 560;
537  cuts[4] = 420;
538  cuts[5] = -8500;
539  cuts[6] = 0.9;
540  }
541  else if (theD0pt > 5){
542  cuts[0] = 800;
543  cuts[1] = 1.0;
544  cuts[2] = 1.2;
545  cuts[3] = 700;
546  cuts[4] = 700;
547  cuts[5] = 10000;
548  cuts[6] = 0.9;
549  }
550  if (dca < cuts[0]
551  && TMath::Abs(cosThetaStar) < cuts[1]
552  && pTpi > cuts[2]
553  && pTK > cuts[2]
554  && TMath::Abs(d0pi) < cuts[3]
555  && TMath::Abs(d0K) < cuts[4]
556  && d0xd0 < cuts[5]
557  && cosPointingAngle > cuts[6]
558  ){
559 
560  AliDebug(2,"Particle passed D* selection cuts");
561  fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;
562 
563  if(!fAcceptanceUnf){ // unfolding
564 
565  Double_t fill[4]; //fill response matrix
566 
567  // dimensions 0&1 : pt,eta (Rec)
568  fill[0] = pt ;
569  fill[1] = rapidity;
570  // dimensions 2&3 : pt,eta (MC)
571  fill[2] = mcVtxHFDStar->Pt();
572  fill[3] = mcVtxHFDStar->Y();
573 
574  fCorrelation->Fill(fill);
575  }
576  }
577  }
578  }
579  }
580  if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
581  } // end loop on D*->Kpipi
582 
583  fHistEventsProcessed->Fill(0);
584 
585  PostData(1,fHistEventsProcessed) ;
586  PostData(2,fCFManager->GetParticleContainer()) ;
587  PostData(3,fCorrelation) ;
588 }
589 
590 
591 //___________________________________________________________________________
593 {
594  // The Terminate()
595  AliAnalysisTaskSE::Terminate();
596 
597  // draw correlation matrix
598  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
599  if(!cont) {
600  printf("CONTAINER NOT FOUND\n");
601  return;
602  }
603 }
604 
605 //___________________________________________________________________________
607 
608  //
609  // useroutput
610  //
611 
612  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
613 
614  //slot #1
615  OpenFile(1);
616  fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
617 }
618 
619 //________________________________________________________________________________
620 Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
621 
622  //
623  // fill the D* and D0 MC container
624  //
625 
626  Bool_t isDStar = kFALSE;
627 
628  if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
629 
630  // getting the daughters
631  Int_t daughter0 = mcPart->GetDaughter(0);
632  Int_t daughter1 = mcPart->GetDaughter(1);
633 
634  AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
635  if (daughter0 == 0 || daughter1 == 0) {
636  AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
637  return isDStar;
638  }
639 
640  if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
641  AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
642  return isDStar;
643  }
644 
645  AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
646  AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
647  if (!mcPartDaughter0 || !mcPartDaughter1) {
648  AliWarning("D*: At least one Daughter Particle not found in tree, skipping");
649  return isDStar;
650  }
651 
652  if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
653  TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
654  !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
655  TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
656  AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
657  return isDStar;
658  }
659 
660  Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
661  Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
662 
663  // getting vertex from daughters
664  mcPartDaughter0->XvYvZv(vtx2daughter0);
665  mcPartDaughter1->XvYvZv(vtx2daughter1);
666 
667  // check if the secondary vertex is the same for both
668  if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
669  AliError("The D* daughters have different secondary vertex, skipping the track");
670  return isDStar;
671  }
672 
673  AliAODMCParticle* neutralDaugh = mcPartDaughter0;
674 
675  Double_t VectorD0[2] ={0.,0.};
676 
677  if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
678  AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
679  return isDStar;
680  }
681  // get the pT of the daughters
682 
683  Double_t pTpi = 0.;
684  Double_t pTD0 = 0.;
685 
686  if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
687  pTpi = mcPartDaughter0->Pt();
688  pTD0 = mcPartDaughter1->Pt();
689  }
690  else {
691  pTpi = mcPartDaughter1->Pt();
692  pTD0 = mcPartDaughter0->Pt();
693  }
694 
695  vectorMC[0] = mcPart->Pt();
696  vectorMC[1] = mcPart->Y() ;
697  vectorMC[2] = 0;
698  vectorMC[3] = pTpi ;
699  vectorMC[4] = pTD0 ;
700  vectorMC[5] = 0;
701  vectorMC[6] = mcPart->Phi() ;
702  vectorMC[7] = VectorD0[0] ;
703  vectorMC[8] = VectorD0[1] ;
704 
705  isDStar = kTRUE;
706 
707  return isDStar;
708 }
709 //________________________________________________________________________________________________
710 
711 Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
712 
713  //
714  // chack wether D0 is decaing into kpi
715  //
716 
717  Bool_t isHadronic = kFALSE;
718 
719  Int_t daughterD00 = neutralDaugh->GetDaughter(0);
720  Int_t daughterD01 = neutralDaugh->GetDaughter(1);
721 
722  AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
723  if (daughterD00 == 0 || daughterD01 == 0) {
724  AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
725  return isHadronic;
726  }
727 
728  if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
729  AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
730  return isHadronic;
731  }
732 
733  AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
734  AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
735  if (!mcPartDaughterD00 || !mcPartDaughterD01) {
736  AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
737  return isHadronic;
738  }
739 
740  if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
741  TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
742  !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
743  TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
744  AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
745  return isHadronic;
746  }
747 
748  Double_t pTD0pi = 0;
749  Double_t pTD0K = 0;
750 
751 
752  if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
753  pTD0pi = mcPartDaughterD00->Pt();
754  pTD0K = mcPartDaughterD01->Pt();
755  }
756  else {
757  pTD0pi = mcPartDaughterD01->Pt();
758  pTD0K = mcPartDaughterD00->Pt();
759  }
760 
761  isHadronic = kTRUE;
762 
763  VectorD0[0] = pTD0pi;
764  VectorD0[1] = pTD0K;
765 
766  return isHadronic;
767 
768 }
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
Double_t YDstar() const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
char Char_t
Definition: External.C:18
Int_t fCountRecoDStarSel
response matrix for unfolding
TCanvas * c
Definition: TestFitELoss.C:172
Bool_t fAcceptanceUnf
min n. of ITS clusters for RecoDecay soft pion
Bool_t GetDStarMCParticle(AliAODMCParticle *mcPart, TClonesArray *mcArray, Double_t *vectorMC) const
AliCFTaskForDStarAnalysis & operator=(const AliCFTaskForDStarAnalysis &c)
Bool_t EvaluateIfD0toKpi(AliAODMCParticle *neutralDaugh, TClonesArray *mcArray, Double_t *VectorD0) const
int Int_t
Definition: External.C:63
Int_t fMinITSClustersSoft
min n. of ITS clusters for RecoDecay
Definition: External.C:204
Double_t CosThetaStarD0bar() const
angle of K
Int_t fEvents
Reco particle found that satisfy cuts in D* selection.
AliAODTrack * GetBachelor() const
AliAODVertex * GetOwnPrimaryVtx() const
rapidity
Definition: HFPtSpectrum.C:47
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
void UserCreateOutputObjects()
ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects.
TH1I * fHistEventsProcessed
pointer to the CF manager
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
AliAODRecoDecayHF2Prong * Get2Prong() const
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65