AliPhysics  7273240 (7273240)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFlavourJetCorrelations.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 //
17 // Analysis Taks for reconstructed particle correlation
18 // (first implementation done for D mesons) with jets
19 // (use the so called Emcal framework)
20 //
21 //-----------------------------------------------------------------------
22 // Authors:
23 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
24 // S. Antônio (University of São Paulo) antonio.silva@cern.ch
25 // A. Grelli (Utrecht University) a.grelli@uu.nl
26 // X. Zhang (LBNL) XMZhang@lbl.gov
27 //-----------------------------------------------------------------------
28 
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
31 #include "TROOT.h"
32 #include <THnSparse.h>
33 #include <TSystem.h>
34 #include <TObjectTable.h>
35 #include "AliMultSelection.h"
36 
38 #include "AliAODHandler.h"
39 #include "AliAnalysisManager.h"
40 #include "AliEmcalJet.h"
41 #include "AliJetContainer.h"
42 #include "AliAODRecoDecay.h"
43 #include "AliAODRecoCascadeHF.h"
45 #include "AliAODMCParticle.h"
46 #include "AliRDHFCutsD0toKpi.h"
48 #include "AliParticleContainer.h"
49 #include "AliEmcalParticle.h"
50 
52 
53 
54 //_______________________________________________________________________________
55 
58 fUseMCInfo(kTRUE),
59 fUseReco(kTRUE),
60 fCandidateType(),
61 fCorrelationMethod(),
62 fPDGmother(),
63 fNProngs(),
64 fPDGdaughters(),
65 fBranchName(),
66 fCuts(0),
67 fMinMass(),
68 fMaxMass(),
69 fCandidateArray(0),
70 fSideBandArray(0),
71 fAnalyseDBkg(kFALSE),
72 fNAxesBigSparse(9),
73 fUseCandArray(kFALSE),
74 fUseSBArray(kFALSE),
75 fhstat(),
76 fhCentDjet(),
77 fhPtJetTrks(),
78 fhPhiJetTrks(),
79 fhEtaJetTrks(),
80 fhPtJet(),
81 fhPhiJet(),
82 fhEtaJet(),
83 fhInvMassptD(),
84 fhDiffSideBand(),
85 fhInvMassptDbg(),
86 fhPtPion(),
87 fhsDphiz(),
88 fResponseMatrix()
89 
90 {
91  //
92  // Default ctor
93 }
94 
95 //_______________________________________________________________________________
96 
98 AliAnalysisTaskEmcalJet(name,kTRUE),
99 fUseMCInfo(kTRUE),
100 fUseReco(kTRUE),
101 fCandidateType(),
102 fCorrelationMethod(),
103 fPDGmother(),
104 fNProngs(),
105 fPDGdaughters(),
106 fBranchName(),
107 fCuts(0),
108 fMinMass(),
109 fMaxMass(),
110 fCandidateArray(0),
111 fSideBandArray(0),
112 fAnalyseDBkg(kFALSE),
113 fNAxesBigSparse(9),
114 fUseCandArray(kFALSE),
115 fUseSBArray(kFALSE),
116 fhstat(),
117 fhCentDjet(),
118 fhPtJetTrks(),
119 fhPhiJetTrks(),
120 fhEtaJetTrks(),
121 fhPtJet(),
122 fhPhiJet(),
123 fhEtaJet(),
124 fhInvMassptD(),
125 fhDiffSideBand(),
126 fhInvMassptDbg(),
127 fhPtPion(),
128 fhsDphiz(),
129 fResponseMatrix()
130 {
131  //
132  // Constructor. Initialization of Inputs and Outputs
133  //
134 
135  Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
136  fCuts=cuts;
137  fCandidateType=candtype;
138  const Int_t nptbins=fCuts->GetNPtBins();
139  Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};
140 
141  switch(fCandidateType){
142  case 0:
143  fPDGmother=421;
144  fNProngs=2;
145  fPDGdaughters[0]=211;//pi
146  fPDGdaughters[1]=321;//K
147  fPDGdaughters[2]=0; //empty
148  fPDGdaughters[3]=0; //empty
149  fBranchName="D0toKpi";
150  break;
151  case 1:
152  fPDGmother=413;
153  fNProngs=3;
154  fPDGdaughters[1]=211;//pi soft
155  fPDGdaughters[0]=421;//D0
156  fPDGdaughters[2]=211;//pi fromD0
157  fPDGdaughters[3]=321; // K from D0
158  fBranchName="Dstar";
159 
160  if(nptbins<=13){
161  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
162  } else {
163  AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
164  }
165  break;
166  default:
167  printf("%d not accepted!!\n",fCandidateType);
168  break;
169  }
170 
173  if(fUseCandArray) DefineInput(1, TClonesArray::Class());
174  if(fUseSBArray) DefineInput(2, TClonesArray::Class());
175 
176  DefineOutput(1,TList::Class()); // histos
177  DefineOutput(2,AliRDHFCuts::Class()); // my cuts
178 
179 }
180 
181 //_______________________________________________________________________________
182 
184  //
185  // destructor
186  //
187 
188  Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
189 
190  delete fCuts;
191 
192 }
193 
194 //_______________________________________________________________________________
195 
197  //
198  // Initialization
199  //
200 
201  if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
202 
203  switch(fCandidateType){
204  case 0:
205  {
206  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
207  copyfCuts->SetName("AnalysisCutsDzero");
208  // Post the data
209  PostData(2,copyfCuts);
210  }
211  break;
212  case 1:
213  {
214  AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
215  copyfCuts->SetName("AnalysisCutsDStar");
216  // Post the cuts
217  PostData(2,copyfCuts);
218  }
219  break;
220  default:
221  return;
222  }
223 
224  return;
225 }
226 
227 //_______________________________________________________________________________
228 
230  // output
231  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
233 
234  // define histograms
235  // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
237  PostData(1,fOutput);
238 
239  return;
240 }
241 
242 //_______________________________________________________________________________
243 
245 {
246  // user exec from AliAnalysisTaskEmcal is used
247 
248  // Load the event
249  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
250 
251  TClonesArray *arrayDStartoD0pi=0;
252 
253  if (!aodEvent && AODEvent() && IsStandardAOD()) {
254 
255  // In case there is an AOD handler writing a standard AOD, use the AOD
256  // event in memory rather than the input (ESD) event.
257  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
258 
259  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
260  // have to taken from the AOD event hold by the AliAODExtension
261  AliAODHandler* aodHandler = (AliAODHandler*)
262  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
263  if(aodHandler->GetExtensions()) {
264  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
265  AliAODEvent *aodFromExt = ext->GetAOD();
266  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
267  }
268  } else if(aodEvent){
269  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
270  }
271 
272  if (!arrayDStartoD0pi) {
273  AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
274  // return;
275  } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
276 
277  TClonesArray* mcArray = 0x0;
278  if (fUseMCInfo) { //not used at the moment,uncomment return if you use
279  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
280  if (!mcArray) {
281  printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
282  //return kFALSE;
283  }
284  }
285 
286  //D meson candidates. Also background if is MC
287 
288  if(fUseCandArray)
289  {
290  fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
291  if (!fCandidateArray) return kFALSE;
292  for(Int_t icand=0; icand<fCandidateArray->GetEntriesFast(); icand++)
293  {
294  fhstat->Fill(2);
295  }
296  }
297  if(fUseSBArray)
298  {
299  fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
300  if (!fSideBandArray) return kFALSE;
301  }
302 
303  fhstat->Fill(0);
304 
305  // fix for temporary bug in ESDfilter
306  // the AODs with null vertex pointer didn't pass the PhysSel
307  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
308 
309  //Event selection
310  Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
311  TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
312  if(!iseventselected) return kFALSE;
313 
314  fhstat->Fill(1);
315  Float_t lPercentile = 300;
316  AliMultSelection *MultSelection = 0x0;
317  MultSelection = (AliMultSelection*)aodEvent->FindListObject("MultSelection");
318  if(!MultSelection) {
319  //If you get this warning (and lPercentiles 300) please check that the AliMultSelectionTask actually ran (before your task)
320  AliWarning("AliMultSelection object not found!");
321  }else{
322  lPercentile = MultSelection->GetMultiplicityPercentile("V0M");
323  }
324 
325  fhCentDjet->Fill(lPercentile);
326 
327  AliJetContainer* JetCont = GetJetContainer(0);
328  if(!JetCont) return kFALSE;
329  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
330 
331  AliJetContainer* JetContSB = 0;
332  AliParticleContainer *ParticlesContSB = 0;
333  if(fAnalyseDBkg)
334  {
335  JetContSB = GetJetContainer(1);
336  if(!JetContSB) return kFALSE;
337  ParticlesContSB = JetContSB->GetParticleContainer();
338  }
339 
340  //Distribution of all particles in the event
341  Int_t ntrarr=ParticlesCont->GetNParticles();
342  for(Int_t i=0;i<ntrarr;i++)
343  {
344  AliVParticle* jetTrk= ParticlesCont->GetParticle(i);
345  if (!jetTrk) continue;
346  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
347  if (emcpart) jetTrk = emcpart->GetTrack();
348  fhPtJetTrks->Fill(jetTrk->Pt());
349  fhPhiJetTrks->Fill(jetTrk->Phi());
350  fhEtaJetTrks->Fill(jetTrk->Eta());
351  }
352 
353  JetCont->ResetCurrentID();
354  AliEmcalJet* jet=0;
355 
356  while ((jet = JetCont->GetNextJet()))
357  {
358  UInt_t rejectionReason = 0;
359  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
360  if(!OKjet) {
361  fhstat->Fill(5);
362  continue;
363  }
364 
365 
366  Double_t JetPtCorr = jet->Pt() - jet->Area()*JetCont->GetRhoVal(); //background subtraction
367  fhstat->Fill(3); //Jet accepted
368  fhPhiJet->Fill(jet->Phi());
369  fhEtaJet->Fill(jet->Eta());
370  fhPtJet->Fill(JetPtCorr);
371  }
372 
373  if(ParticlesCont->GetNParticles()>0) fhstat->Fill(2);
374 
376  {
377  if(ParticlesCont->GetNParticles()>0) ConstituentCorrelationMethod(kFALSE,aodEvent);
378  if(fAnalyseDBkg==kTRUE && ParticlesContSB->GetNParticles()>0) ConstituentCorrelationMethod(kTRUE,aodEvent);
379  }
380 
381  else if(fCorrelationMethod==kAngular)
382  {
383  if(fCandidateArray->GetEntriesFast()>0) AngularCorrelationMethod(kFALSE,aodEvent);
384  if(fAnalyseDBkg==kTRUE && fSideBandArray->GetEntriesFast()>0) AngularCorrelationMethod(kTRUE,aodEvent);
385  }
387  {
389  }
390 
391  PostData(1,fOutput);
392  return kTRUE;
393 }
395 {
396 
397  AliJetContainer* JetCont = 0;
398 
399  if(!IsBkg) JetCont = GetJetContainer(0);
400  else JetCont = GetJetContainer(1);
401 
402  Double_t rho = 0;
403  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
404 
405  AliEmcalJet *jet;
406 
407  GetHFJet(jet,IsBkg);
408 
409  if(jet) FillDJetHistograms(jet,rho,IsBkg,aodEvent);
410 
411 }
413 {
414  AliJetContainer* JetCont = GetJetContainer(0);
415 
416  Int_t ncand = 0;
417  if(!IsBkg) ncand = fCandidateArray->GetEntriesFast();
418  else ncand = fSideBandArray->GetEntriesFast();
419 
420  Double_t rho = 0;
421  if(!JetCont->GetRhoName().IsNull()) rho = JetCont->GetRhoVal();
422 
423  for(Int_t icand = 0; icand<ncand; icand++)
424  {
425  AliVParticle* charm=0x0;
426  if(!IsBkg) charm = (AliVParticle*)fCandidateArray->At(icand);
427  else charm = (AliVParticle*)fSideBandArray->At(icand);
428  if(!charm) continue;
429 
430  Int_t JetTag = AliEmcalJet::kD0;
432  //loop over jets
433  JetCont->ResetCurrentID();
434  AliEmcalJet* jet=0;
435  while ((jet = JetCont->GetNextJet()))
436  {
437  UInt_t rejectionReason = 0;
438  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
439  if(!OKjet) continue;
440 
441  if(DeltaR(jet,charm,rho)<JetCont->GetJetRadius() && CheckDeltaR(jet,charm)<JetCont->GetJetRadius())
442  {
443  if(!IsBkg) jet->AddFlavourTag(JetTag);
444  jet->AddFlavourTrack(charm);
445  FillDJetHistograms(jet,rho,IsBkg,aodEvent);
446  }
447  }
448  }
449 
450 }
451 
453 {
454  AliJetContainer* JetContRec = GetJetContainer(0);
455  //AliJetContainer* JetContGen = GetJetContainer(1);
456 
457  Double_t rho = JetContRec->GetRhoVal();
458 
459  AliEmcalJet *jetRec = 0;
460  AliEmcalJet *jetGen = 0;
461 
462  GetHFJet(jetRec,kFALSE);
463  GetHFJet(jetGen,kTRUE);
464  //GetHFJetParticleLevel(jetGen,kTRUE);
465 
466  if(!jetRec) AliDebug(2, "No Reconstructed Level Jet Found!");
467  else if(!jetGen) AliDebug(2, "No Generated Level Jet Found!");
468  else
469  {
470  AliVParticle *Drec = jetRec->GetFlavourTrack();
471  AliVParticle *Dgen = jetGen->GetFlavourTrack();
472  Double_t zRec = Z(Drec,jetRec,rho);
473  Double_t zGen = Z(Dgen,jetGen,0);
474  Double_t JetPtRec = jetRec->Pt() - rho*jetRec->Area();
475  Double_t fillRM[6] = {zRec,JetPtRec,Drec->Pt(),zGen,jetGen->Pt(),Dgen->Pt()};
476  fResponseMatrix->Fill(fillRM,1.);
477  }
478 }
479 
480 //_______________________________________________________________________________
481 
483 {
484 
485  AliVParticle *Dmeson = jet->GetFlavourTrack(0);
486  Double_t JetPtCorr = jet->Pt() - rho*jet->Area();
487  Double_t z = 0;
488  if(rho>0) z = Z(Dmeson,jet,rho);
489  else z = Z(Dmeson,jet);
490 
491  Bool_t bDInEMCalAcc=InEMCalAcceptance(Dmeson);
492  Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
493 
495  {
496  AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)Dmeson;
497  FillHistogramsD0JetCorr(dzero,z,Dmeson->Pt(),JetPtCorr,jet->Phi(),IsBkg,bDInEMCalAcc,bJetInEMCalAcc,aodEvent);
498  }
499  else if(fCandidateType==kDstartoKpipi)
500  {
502  FillHistogramsDstarJetCorr(dstar,z,Dmeson->Pt(),JetPtCorr,jet->Phi(),IsBkg,bDInEMCalAcc,bJetInEMCalAcc);
503  }
504 
505 
506 }
507 
509 {
510  AliJetContainer* JetCont = 0;
511 
512  if(!IsBkg) JetCont = GetJetContainer(0);
513  else JetCont = GetJetContainer(1);
514 
515  AliParticleContainer *ParticlesCont = JetCont->GetParticleContainer();
516 
517  Bool_t JetIsHF = kFALSE;
518 
519  JetCont->ResetCurrentID();
520  while ((jet = JetCont->GetNextJet()))
521  {
522  UInt_t rejectionReason = 0;
523  Bool_t OKjet = JetCont->AcceptJet(jet, rejectionReason);
524  if(!OKjet) continue;
525 
526  Int_t JetTag = AliEmcalJet::kD0;
527  TString recoDecayClassName("AliAODRecoDecayHF2Prong");
529  {
530  JetTag = AliEmcalJet::kDStar;
531  recoDecayClassName = "AliAODRecoCascadeHF";
532  }
533  //loop on jet particles
534  Int_t ntrjet= jet->GetNumberOfTracks();
535  for(Int_t itrk=0;itrk<ntrjet;itrk++)
536  {
537  AliVParticle* jetTrk=jet->TrackAt(itrk,ParticlesCont->GetArray());
538  if(!jetTrk) continue;
539  AliEmcalParticle* emcpart = dynamic_cast<AliEmcalParticle*>(jetTrk);
540  if(emcpart) jetTrk = emcpart->GetTrack();
541  if(strcmp(jetTrk->ClassName(),recoDecayClassName)==0)
542  {
543  JetIsHF = kTRUE;
544  if(!IsBkg) jet->AddFlavourTag(JetTag);
545  jet->AddFlavourTrack(jetTrk);
546  break;
547  }
548  } //end loop on jet tracks
549  if(JetIsHF) break;
550  } //end of jet loop
551 
552  if(!JetIsHF) jet = 0;
553 
554 }
555 
556 //_______________________________________________________________________________
557 
559 {
560  // The Terminate() function is the last function to be called during
561  // a query. It always runs on the client, it can be used to present
562  // the results graphically or save the results to file.
563 
564  Info("Terminate"," terminate");
565  AliAnalysisTaskSE::Terminate();
566 
567  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
568  if (!fOutput) {
569  printf("ERROR: fOutput not available\n");
570  return;
571  }
572 }
573 
574 //_______________________________________________________________________________
575 
577  Float_t mass=0;
578  Int_t abspdg=TMath::Abs(pdg);
579 
580  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
581  // compute the Delta mass for the D*
583  Float_t mass1=0;
584  mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
585  mass = mass-mass1;
586  }
587 
588  fMinMass = mass-range;
589  fMaxMass = mass+range;
590 
591  AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
592  if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
593 }
594 
595 //_______________________________________________________________________________
596 
598  if(uplimit>lowlimit)
599  {
600  fMinMass = lowlimit;
601  fMaxMass = uplimit;
602  }
603  else{
604  printf("Error! Lower limit larger than upper limit!\n");
605  fMinMass = uplimit - uplimit*0.2;
606  fMaxMass = uplimit;
607  }
608 }
609 
610 //_______________________________________________________________________________
611 
613  if(nptbins>30) {
614  AliInfo("Maximum number of bins allowed is 30!");
615  return kFALSE;
616  }
617  if(!width) return kFALSE;
618  for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
619  return kTRUE;
620 }
621 
622 //_______________________________________________________________________________
623 
625 
626  Double_t p[3],pj[3];
627  Bool_t okpp=part->PxPyPz(p);
628  Bool_t okpj=jet->PxPyPz(pj);
629 
630  if(!okpp || !okpj)
631  {
632  printf("Problems getting momenta\n");
633  return -999;
634  }
635 
636  //Background Subtraction
637  //It corrects the each component of the jet momentum for Z calculation
638 
639  pj[0] = jet->Px() - jet->Area()*(rho*TMath::Cos(jet->AreaPhi()));
640  pj[1] = jet->Py() - jet->Area()*(rho*TMath::Sin(jet->AreaPhi()));
641  pj[2] = jet->Pz() - jet->Area()*(rho*TMath::SinH(jet->AreaEta()));
642 
643  return Z(p,pj);
644 }
645 //_______________________________________________________________________________
646 
648 
649  Double_t p[3],pj[3];
650  Bool_t okpp=part->PxPyPz(p);
651  Bool_t okpj=jet->PxPyPz(pj);
652 
653  if(!okpp || !okpj)
654  {
655  printf("Problems getting momenta\n");
656  return -999;
657  }
658 
659  return Z(p,pj);
660 }
661 //_______________________________________________________________________________
663 
664  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
665  Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
666  return z;
667 }
668 
669 
670 //_______________________________________________________________________________
672 
673  Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
674  Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
675  return z;
676 }
677 
678 //_______________________________________________________________________________
679 
681 
682  // Statistics
683  Int_t nbins=8;
684  if(fUseMCInfo) nbins+=2;
685 
686  fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
687  fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
688  fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
689  fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
690  fhstat->GetXaxis()->SetBinLabel(4,"N jets");
691  fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
692  fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
693  fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
694  fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
695  if(fUseMCInfo) {
696  fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
697  fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
698  fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
699  fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
700  }
701  fhstat->SetNdivisions(1);
702  fOutput->Add(fhstat);
703 
704  fhCentDjet=new TH1F("hCentDjet","Centrality",100,0,100);
705  fOutput->Add(fhCentDjet);
706 
707  const Int_t nbinsmass=300;
708  const Int_t nbinspttrack=500;
709  Int_t nbinsptjet=200;
710  const Int_t nbinsptD=100;
711  const Int_t nbinsphi=200;
712  const Int_t nbinseta=100;
713 
714  //binning for THnSparse
715  const Int_t nbinsSpsmass=60;
716  const Int_t nbinsSpsptjet=200;
717  const Int_t nbinsSpsptD=50;
718  const Int_t nbinsSpsz=144;
719 
720  const Float_t pttracklims[2]={0.,200.};
721  Float_t ptjetlims[2]={-50.,150.};
722  const Float_t ptDlims[2]={0.,50.};
723  const Float_t zlims[2]={0.,1.2};
724  const Float_t philims[2]={0.,6.3};
725  const Float_t etalims[2]={-1.5,1.5};
726 
727  // jet related fistograms
728 
729  fhPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
730  fhPhiJetTrks->Sumw2();
731  fhEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
732  fhEtaJetTrks->Sumw2();
733  fhPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinspttrack,pttracklims[0],pttracklims[1]);
734  fhPtJetTrks->Sumw2();
735 
736  fhPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);
737  fhPhiJet->Sumw2();
738  fhEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
739  fhEtaJet->Sumw2();
740  fhPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
741  fhPtJet->Sumw2();
742 
743  fOutput->Add(fhPhiJetTrks);
744  fOutput->Add(fhEtaJetTrks);
745  fOutput->Add(fhPtJetTrks);
746  fOutput->Add(fhPhiJet);
747  fOutput->Add(fhEtaJet);
748  fOutput->Add(fhPtJet);
749 
751  {
752  if(fAnalyseDBkg){
753  fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
754  fhDiffSideBand->SetStats(kTRUE);
755  fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
756  fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
757  fhDiffSideBand->Sumw2();
758  fOutput->Add(fhDiffSideBand);
759  }
760 
761  fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
762  fhPtPion->SetStats(kTRUE);
763  fhPtPion->GetXaxis()->SetTitle("GeV/c");
764  fhPtPion->GetYaxis()->SetTitle("Entries");
765  fhPtPion->Sumw2();
766  fOutput->Add(fhPtPion);
767 
768  }
769 
770 
771  fhInvMassptD = new TH2F("hInvMassptD","D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
772  fhInvMassptD->SetStats(kTRUE);
773  fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
774  fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
775  fhInvMassptD->Sumw2();
776 
777  fOutput->Add(fhInvMassptD);
778 
779  if(fUseMCInfo){
780  fhInvMassptDbg = new TH2F("hInvMassptDbg","Background D (Delta R < Rjet) invariant mass distribution p_{T}^{j} > threshold",nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
781  fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
782  fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
783  fhInvMassptDbg->Sumw2();
784  fOutput->Add(fhInvMassptDbg);
785 
786  }
787 
788  fhsDphiz=0x0; //definition below according to the switches
789 
790  if(fUseMCInfo){
791  AliInfo("Creating a 7 axes container (MB background candidates)");
792  const Int_t nAxis=7;
793  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
794  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
795  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
796  fNAxesBigSparse=nAxis;
797  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
798 
799  } else{
800  AliInfo("Creating a 5 axes container");
801  const Int_t nAxis=5;
802  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,63};
803  const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,0};
804  const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass,6.3};
805  fNAxesBigSparse=nAxis;
806 
807  fhsDphiz=new THnSparseF("hsDphiz","Z, p_{T}^{jet}, p_{T}^{D}, mass., Jet phi", nAxis, nbinsSparse, minSparse, maxSparse);
808  }
809 
810  if(!fhsDphiz) AliFatal("No THnSparse created");
811  fhsDphiz->Sumw2();
812 
813  fOutput->Add(fhsDphiz);
814 
816  {
817  const Int_t RMnbins[6] = {nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsz,nbinsSpsptjet,nbinsSpsptD};
818  const Double_t RMmin[6]={zlims[0],ptjetlims[0],ptDlims[0],zlims[0],ptjetlims[0],ptDlims[0]};
819  const Double_t RMmax[6]={zlims[1],ptjetlims[1],ptDlims[1],zlims[1],ptjetlims[1],ptDlims[1]};
820  fResponseMatrix = new THnSparseF("ResponseMatrix","z, jet pt, D pt: Rec and Gen",6,RMnbins,RMmin,RMmax);
821  fResponseMatrix->Sumw2();
822  fOutput->Add(fResponseMatrix);
823  }
824 
825  PostData(1, fOutput);
826 
827  return kTRUE;
828 }
829 
830 //_______________________________________________________________________________
831 
832 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent){
833 
834 
835  Double_t masses[2]={0.,0.};
836  Int_t pdgdaughtersD0[2]={211,321};//pi,K
837  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
838 
839  masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
840  masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
841 
842  Double_t *point=0x0;
843 
844  if(!fAnalyseDBkg)
845  {
846  point=new Double_t[5];
847  point[0]=z;
848  point[1]=ptj;
849  point[2]=ptD;
850  point[3]=masses[0];
851  point[4]=phij;
852 
853  }
854  else
855  {
856  point=new Double_t[7];
857  point[0]=z;
858  point[1]=ptj;
859  point[2]=ptD;
860  point[3]=masses[0];
861  point[4]=static_cast<Double_t>(IsBkg ? 1 : 0);
862  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
863  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
864 
865  }
866  Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
867  if(isselected==1 || isselected==3)
868  {
869 
870  fhInvMassptD->Fill(masses[0],ptD);
871  fhsDphiz->Fill(point,1.);
872 
873  }
874  if(isselected>=2)
875  {
876  fhInvMassptD->Fill(masses[1],ptD);
877  point[3]=masses[1];
878  fhsDphiz->Fill(point,1.);
879 
880  }
881  delete[] point;
882 }
883 
884 //_______________________________________________________________________________
885 
887 
888  AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
889  Double_t deltamass= dstar->DeltaInvMass();
890 
891 
892  fhPtPion->Fill(softpi->Pt());
893 
894  Double_t *point=0x0;
895  if(!fAnalyseDBkg)
896  {
897  point=new Double_t[5];
898  point[0]=z;
899  point[1]=ptj;
900  point[2]=ptD;
901  point[3]=deltamass;
902  point[4]=phij;
903  }
904  else
905  {
906  point=new Double_t[7];
907  point[0]=z;
908  point[1]=ptj;
909  point[2]=ptD;
910  point[3]=deltamass;
911  point[4]=static_cast<Double_t>(IsBkg ? 1 : 0);
912  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
913  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
914  }
915 
916  if(!point){
917  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
918  return;
919  }
920 
921 
922  fhInvMassptD->Fill(deltamass,ptD);
923  fhsDphiz->Fill(point,1.);
924  delete[] point;
925 }
926 
927 //_______________________________________________________________________________
928 
930 
931  Double_t pdgmass=0;
932  if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
933  if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
934 
935  Double_t *point=0x0;
936 
937  if(fNAxesBigSparse==5){
938  point=new Double_t[5];
939  point[0]=z;
940  point[1]=ptjet;
941  point[2]=ptD;
942  point[3]=pdgmass;
943  point[4]=phij;
944  }
945  if(fNAxesBigSparse==7){
946  point=new Double_t[7];
947  point[0]=z;
948  point[1]=ptjet;
949  point[2]=ptD;
950  point[3]=pdgmass;
951  point[4]=1;
952  point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
953  point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
954  }
955 
956  if(!point){
957  AliError(Form("Numer of THnSparse entries %d not valid", fNAxesBigSparse));
958  return;
959  }
960 
961 
962  point[3]=pdgmass;
963  fhsDphiz->Fill(point,1.);
964 
965  delete[] point;
966 }
967 
968 //_______________________________________________________________________________
969 
971  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
972  //It recalculates the eta-phi values if it was asked for background subtraction of the jets
973  if(!p1 || !p2) return -1;
974 
975  Double_t phi1=p1->Phi(),eta1=p1->Eta();
976 
977  if (rho>0)
978  {
979  Double_t pj[3];
980  Bool_t okpj=p1->PxPyPz(pj);
981  if(!okpj){
982  printf("Problems getting momenta\n");
983  return -999;
984  }
985  pj[0] = p1->Px() - p1->Area()*(rho*TMath::Cos(p1->AreaPhi()));
986  pj[1] = p1->Py() - p1->Area()*(rho*TMath::Sin(p1->AreaPhi()));
987  pj[2] = p1->Pz() - p1->Area()*(rho*TMath::SinH(p1->AreaEta()));
988  //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
989  //[0,pi] if py > 0 and
990  //[pi,2pi] if py < 0
991  phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
992  if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
993  eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
994  }
995 
996  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
997 
998  Double_t dPhi=TMath::Abs(phi1-phi2);
999  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1000 
1001 
1002  Double_t dEta=eta1-eta2;
1003  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1004  return deltaR;
1005 
1006 }
1007 //_______________________________________________________________________________
1009  //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1010  if(!p1 || !p2) return -1;
1011 
1012  Double_t phi1=p1->Phi(),eta1=p1->Eta();
1013 
1014  Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1015 
1016  Double_t dPhi=TMath::Abs(phi1-phi2);
1017  if(dPhi>TMath::Pi()) dPhi = (2*TMath::Pi())-dPhi;
1018 
1019 
1020  Double_t dEta=eta1-eta2;
1021  Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1022  return deltaR;
1023 
1024 }
1025 
1026 //_______________________________________________________________________________
1027 
1029 
1030  Double_t ptD=candDstar->Pt();
1031  Int_t bin = fCuts->PtBin(ptD);
1032  if (bin < 0){
1033  // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1034  bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1035  return -1; // out of bounds
1036  }
1037 
1038  Double_t invM=candDstar->InvMassD0();
1039  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1040 
1041  Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1042  Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1043 
1044  if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1045  else return 0;
1046 
1047 }
1048 
1049 //_______________________________________________________________________________
1050 
1052  //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1053 
1054  Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1055  Bool_t binEMCal=kTRUE;
1056  Double_t phi=vpart->Phi(), eta=vpart->Eta();
1057  if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1058  if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1059  return binEMCal;
1060 
1061 
1062 }
Int_t pdg
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:237
Double_t Area() const
Definition: AliEmcalJet.h:117
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
Double_t DeltaInvMass() const
Jet is tagged to contain a D* meson.
Definition: AliEmcalJet.h:72
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:108
Double_t Z(AliVParticle *part, AliEmcalJet *jet, Double_t rho) const
Double_t Py() const
Definition: AliEmcalJet.h:94
Double_t Phi() const
Definition: AliEmcalJet.h:104
Double_t mass
Int_t GetNParticles() const
char Char_t
Definition: External.C:18
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
AliVTrack * GetTrack() const
void AngularCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
AliVParticle * GetFlavourTrack(Int_t i=0) const
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
Double_t Px() const
Definition: AliEmcalJet.h:93
AliParticleContainer * GetParticleContainer() const
Double_t AreaEta() const
Definition: AliEmcalJet.h:119
int Int_t
Definition: External.C:63
void FillDJetHistograms(AliEmcalJet *jet, Double_t rho, Bool_t IsBkg, AliAODEvent *aodEvent)
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void AddFlavourTrack(AliVParticle *hftrack)
virtual AliVParticle * GetParticle(Int_t i=-1) const
AliAODTrack * GetBachelor() const
Float_t CheckDeltaR(AliEmcalJet *p1, AliVParticle *p2) const
ClassImp(AliAnalysisTaskFlavourJetCorrelations) AliAnalysisTaskFlavourJetCorrelations
Double_t AreaPhi() const
Definition: AliEmcalJet.h:120
Double_t Pt() const
Definition: AliEmcalJet.h:96
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:23
AliEmcalJet * GetNextJet()
Bool_t IsEventSelected(AliVEvent *event)
Jet is tagged to contain a D0 meson.
Definition: AliEmcalJet.h:73
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
void FillHistogramsMCGenDJetCorr(Double_t z, Double_t ptD, Double_t ptjet, Double_t phij, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
void ConstituentCorrelationMethod(Bool_t IsBkg, AliAODEvent *aodEvent)
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
void FillHistogramsD0JetCorr(AliAODRecoDecayHF *candidate, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent *aodEvent)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:276
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Bool_t PxPyPz(Double_t p[3]) const
Definition: AliEmcalJet.h:98
Bool_t SetD0WidthForDStar(Int_t nptbins, Float_t *width)
Double_t Pz() const
Definition: AliEmcalJet.h:95
void FillHistogramsDstarJetCorr(AliAODRecoCascadeHF *dstar, Double_t z, Double_t ptD, Double_t ptj, Double_t phij, Bool_t IsBkg, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc)
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:234
const Int_t nbins
bool Bool_t
Definition: External.C:53
Int_t PtBin(Double_t pt) const
Float_t DeltaR(AliEmcalJet *p1, AliVParticle *p2, Double_t rho) const
Int_t nptbins
Container for jet within the EMCAL jet framework.
TClonesArray * fSideBandArray
contains candidates selected by AliRDHFCuts
Bool_t fAnalyseDBkg
contains candidates selected by AliRDHFCuts::IsSelected(kTracks), to be used for side bands (DStar ca...