f28bcf4cfbb38a90635d7a1ba8dcb2badf640e3d
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskFlavourJetCorrelations.cxx
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 // A. Grelli (Utrecht University) a.grelli@uu.nl
25 // X. Zhang (LBNL)  XMZhang@lbl.gov
26 //-----------------------------------------------------------------------
27
28 #include <TDatabasePDG.h>
29 #include <TParticle.h>
30 #include "TROOT.h"
31 #include <TH3F.h>
32 #include <THnSparse.h>
33
34 #include "AliAnalysisTaskFlavourJetCorrelations.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODHandler.h"
37 #include "AliAnalysisManager.h"
38 #include "AliLog.h"
39 #include "AliEmcalJet.h"
40 #include "AliJetContainer.h"
41 #include "AliAODRecoDecay.h"
42 #include "AliAODRecoCascadeHF.h"
43 #include "AliAODRecoDecayHF2Prong.h"
44 #include "AliESDtrack.h"
45 #include "AliAODMCParticle.h"
46 #include "AliPicoTrack.h"
47 #include "AliRDHFCutsD0toKpi.h"
48 #include "AliRDHFCutsDStartoKpipi.h"
49
50 ClassImp(AliAnalysisTaskFlavourJetCorrelations)
51
52
53 //_______________________________________________________________________________
54
55 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
56 AliAnalysisTaskEmcalJet("",kFALSE),
57 fUseMCInfo(kTRUE), 
58 fUseReco(kTRUE),
59 fCandidateType(),
60 fPDGmother(),
61 fNProngs(),
62 fPDGdaughters(),
63 fBranchName(),
64 fmyOutput(0),
65 fCuts(0),
66 fMinMass(),
67 fMaxMass(),  
68 fJetArrName(0),
69 fCandArrName(0),
70 fLeadingJetOnly(kFALSE),
71 fJetRadius(0)
72 {
73    //
74    // Default ctor
75    //
76    //SetMakeGeneralHistograms(kTRUE);
77    
78 }
79
80 //_______________________________________________________________________________
81
82 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
83 AliAnalysisTaskEmcalJet(name,kFALSE),
84 fUseMCInfo(kTRUE),
85 fUseReco(kTRUE),  
86 fCandidateType(),
87 fPDGmother(),
88 fNProngs(),
89 fPDGdaughters(),
90 fBranchName(),
91 fmyOutput(0),
92 fCuts(0),
93 fMinMass(),
94 fMaxMass(),  
95 fJetArrName(0),
96 fCandArrName(0),
97 fLeadingJetOnly(kFALSE),
98 fJetRadius(0)
99 {
100    //
101    // Constructor. Initialization of Inputs and Outputs
102    //
103    
104    Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
105    fCuts=cuts;
106    fCandidateType=candtype;
107    const Int_t nptbins=fCuts->GetNPtBins();
108    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};
109    
110    switch(fCandidateType){
111    case 0:
112       fPDGmother=421;
113       fNProngs=2;
114       fPDGdaughters[0]=211;//pi 
115       fPDGdaughters[1]=321;//K
116       fPDGdaughters[2]=0; //empty
117       fPDGdaughters[3]=0; //empty
118       fBranchName="D0toKpi";
119       fCandArrName="D0";
120       break;
121    case 1: 
122       fPDGmother=413;
123       fNProngs=3;
124       fPDGdaughters[1]=211;//pi soft
125       fPDGdaughters[0]=421;//D0
126       fPDGdaughters[2]=211;//pi fromD0
127       fPDGdaughters[3]=321; // K from D0
128       fBranchName="Dstar";
129       fCandArrName="Dstar";
130       
131       if(nptbins<=13){
132          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
133       } else {
134          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
135       }
136       break;
137    default:
138       printf("%d not accepted!!\n",fCandidateType);
139       break;
140    }
141    
142    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
143    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
144    
145    DefineOutput(1,TList::Class()); // histos
146    DefineOutput(2,AliRDHFCuts::Class()); // my cuts
147    
148 }
149
150 //_______________________________________________________________________________
151
152 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
153    //
154    // destructor
155    //
156    
157    Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");  
158    
159    delete fmyOutput;
160    delete fCuts;
161    
162 }
163
164 //_______________________________________________________________________________
165
166 void AliAnalysisTaskFlavourJetCorrelations::Init(){
167    //
168    // Initialization
169    //
170    
171    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
172    switch(fCandidateType){
173    case 0:
174       {
175          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
176          copyfCuts->SetName("AnalysisCutsDzero");
177          // Post the data
178          PostData(2,copyfCuts);
179       }
180       break;
181    case 1:
182       {
183          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
184          copyfCuts->SetName("AnalysisCutsDStar");
185          // Post the cuts
186          PostData(2,copyfCuts);
187       }
188       break;
189    default:
190       return;
191    }
192    
193    return;
194 }
195
196 //_______________________________________________________________________________
197
198 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { 
199    // output 
200    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
201    fmyOutput = new TList();
202    fmyOutput->SetOwner();
203    fmyOutput->SetName("pippo");
204    // define histograms
205    DefineHistoForAnalysis();
206    PostData(1,fmyOutput);
207    
208    return;
209 }
210
211 //_______________________________________________________________________________
212
213 void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)
214 {
215    // user exec
216    if (!fInitialized){
217       AliAnalysisTaskEmcalJet::ExecOnce();
218    }
219    
220    // Load the event
221    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
222    
223    TClonesArray *arrayDStartoD0pi=0;
224    TClonesArray *trackArr = 0;
225    TClonesArray *candidatesArr = 0;
226    TClonesArray *sbcandArr = 0;
227    
228    if (!aodEvent && AODEvent() && IsStandardAOD()) {
229       
230       // In case there is an AOD handler writing a standard AOD, use the AOD 
231       // event in memory rather than the input (ESD) event.    
232       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
233       
234       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
235       // have to taken from the AOD event hold by the AliAODExtension
236       AliAODHandler* aodHandler = (AliAODHandler*) 
237       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
238       if(aodHandler->GetExtensions()) {
239          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
240          AliAODEvent *aodFromExt = ext->GetAOD();
241          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
242       }
243    } else {
244       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
245    }
246    
247    if (!arrayDStartoD0pi) {
248       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
249       //  return;
250    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
251    
252    TClonesArray* mcArray = 0x0;
253    if (fUseMCInfo) {
254       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
255       if (!mcArray) {
256          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
257          return;
258       }
259    }
260    
261    //retrieve jets
262    trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
263    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
264    //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
265    fJetRadius=GetJetContainer(0)->GetJetRadius();
266    
267    if(!trackArr){
268       AliInfo(Form("Could not find the track array\n"));
269       return;
270    }
271    
272    
273    TString arrname="fCandidateArray";
274    TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
275    candidatesArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));
276    if (!candidatesArr) {
277       Printf("%s array not found",candarrname.Data());
278       InputEvent()->GetList()->ls();
279       return;
280    }
281    //Printf("ncandidates found %d",candidatesArr->GetEntriesFast());
282    
283    TString arrSBname="fSideBandArray";
284    TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
285    sbcandArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));
286    if (fCandidateType==1 && !sbcandArr) {
287       Printf("%s array not found",sbarrname.Data());
288       InputEvent()->GetList()->ls();
289       return;
290    }
291    //Printf("nSBCand or Bkg found %d",sbcandArr->GetEntriesFast());
292    
293    
294    //Histograms
295    TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
296    TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
297    TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
298    TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
299    TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
300    TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
301    TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
302    TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
303    TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
304    TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
305    TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
306    TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
307    TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
308    TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
309    TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
310    TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
311    THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
312        
313    hstat->Fill(0);
314    
315    // fix for temporary bug in ESDfilter 
316    // the AODs with null vertex pointer didn't pass the PhysSel
317    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
318    
319    //Event selection
320    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
321    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
322    if(!iseventselected) return;
323    
324    hstat->Fill(1);
325
326    //retrieve charm candidates selected
327    Int_t candidates = candidatesArr->GetEntriesFast();
328   
329    //trigger on jets
330    
331    Int_t njets=GetJetContainer()->GetNJets();
332    if(njets == 0) {
333       hstat->Fill(6, candidates);
334       hNDPerEvNoJet->Fill(candidates);
335       for(Int_t iD=0;iD<candidates;iD++){
336          AliVParticle* cand=(AliVParticle*)candidatesArr->At(iD);
337          if(!cand) continue;
338          hptDPerEvNoJet->Fill(cand->Pt());
339       
340       }
341       return;
342       
343    }
344    
345    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
346    
347    for(Int_t ic = 0; ic < candidates; ic++) {
348       
349       // D* candidates
350       AliVParticle* charm=0x0;
351       charm=(AliVParticle*)candidatesArr->At(ic);
352       if(!charm) continue;
353       hstat->Fill(2);
354       
355       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
356       
357       if(fCandidateType==kDstartoKpipi) {
358          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
359          Double_t deltamass= dstar->DeltaInvMass();
360          candsparse[3]=deltamass;
361          hnspDstandalone->Fill(candsparse);
362       }
363       if(fCandidateType==kD0toKpi){
364          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
365          Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
366          
367          Double_t masses[2];
368          Int_t pdgdaughtersD0[2]={211,321};//pi,K 
369          Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
370          
371          masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
372          masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
373          if(isselected==1 || isselected==3) {
374             candsparse[3]=masses[0];
375             hnspDstandalone->Fill(candsparse);
376          }
377          if(isselected>=2){
378             candsparse[3]=masses[1];
379             hnspDstandalone->Fill(candsparse);
380             
381          }
382       }
383    }
384
385    // we start with jets
386    Double_t ejet   = 0;
387    Double_t phiJet = 0;
388    Double_t etaJet = 0;
389    Double_t ptjet = 0;
390    Double_t leadingJet =0;
391    
392    Int_t ntrarr=trackArr->GetEntriesFast();
393    hNtrArr->Fill(ntrarr);
394    
395    for(Int_t i=0;i<ntrarr;i++){
396       AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
397       //reusing histograms
398       hPtJetTrks->Fill(jtrack->Pt());
399       hPhiJetTrks->Fill(jtrack->Phi());
400       hEtaJetTrks->Fill(jtrack->Eta());
401       hEjetTrks->Fill(jtrack->E());
402    }
403    
404    
405    //option to use only the leading jet
406    if(fLeadingJetOnly){
407       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {    
408          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
409          ptjet   = jetL->Pt();
410          if(ptjet>leadingJet ) leadingJet = ptjet;
411          
412       }
413    }
414    
415    Int_t cntjet=0;
416    //loop over jets and consider only the leading jet in the event
417    for (Int_t iJets = 0; iJets<njets; iJets++) {    
418       //Printf("Jet N %d",iJets);
419       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
420       //Printf("Corr task Accept Jet");
421       if(!AcceptJet(jet)) {
422          hstat->Fill(5);
423          continue;
424       }
425       
426       //jets variables
427       ejet   = jet->E();
428       phiJet = jet->Phi();
429       etaJet = jet->Eta();
430       ptjet = jet->Pt();
431       
432       // choose the leading jet
433       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
434       //used for normalization
435       hstat->Fill(3);
436       cntjet++;
437       // fill energy, eta and phi of the jet
438       hEjet   ->Fill(ejet);
439       hPhiJet ->Fill(phiJet);
440       hEtaJet ->Fill(etaJet);
441       hPtJet  ->Fill(ptjet);
442       
443       //loop on jet particles
444       Int_t ntrjet=  jet->GetNumberOfTracks();    
445       for(Int_t itrk=0;itrk<ntrjet;itrk++){
446          
447          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);     
448          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
449          
450       }//end loop on jet tracks
451       
452       if(candidates==0){
453          hstat->Fill(7);
454          hPtJetPerEvNoD->Fill(jet->Pt());
455       }
456       
457       //Printf("N candidates %d ", candidates);
458       for(Int_t ic = 0; ic < candidates; ic++) {
459          
460          // D* candidates
461          AliVParticle* charm=0x0;
462          charm=(AliVParticle*)candidatesArr->At(ic);
463          if(!charm) continue;
464          
465          
466          FlagFlavour(charm, jet);
467          if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
468          
469          FillHistogramsRecoJetCorr(charm, jet);
470          
471       }
472       //retrieve side band background candidates for Dstar
473       Int_t nsbcand = 0; if(sbcandArr) nsbcand=sbcandArr->GetEntriesFast();
474       
475       for(Int_t ib=0;ib<nsbcand;ib++){
476          if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
477             AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)sbcandArr->At(ib);
478             if(!sbcand) continue;
479             SideBandBackground(sbcand,jet);
480          }
481          if(fUseMCInfo){
482             AliAODRecoDecayHF* charmbg = 0x0;
483             charmbg=(AliAODRecoDecayHF*)candidatesArr->At(ib);
484             if(!charmbg) continue;
485             MCBackground(charmbg,jet);
486          }
487       }
488    } // end of jet loop
489    
490    hNJetPerEv->Fill(cntjet);
491    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
492    PostData(1,fmyOutput);
493    
494 }
495
496 //_______________________________________________________________________________
497
498 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
499 {    
500    // The Terminate() function is the last function to be called during
501    // a query. It always runs on the client, it can be used to present
502    // the results graphically or save the results to file.
503    
504    Info("Terminate"," terminate");
505    AliAnalysisTaskSE::Terminate();
506    
507    fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
508    if (!fmyOutput) {     
509       printf("ERROR: fmyOutput not available\n");
510       return;
511    }
512 }
513
514 //_______________________________________________________________________________
515
516 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
517    Float_t mass=0;
518    Int_t abspdg=TMath::Abs(pdg);
519    
520    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
521    // compute the Delta mass for the D*
522    if(fCandidateType==kDstartoKpipi){
523       Float_t mass1=0;
524       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
525       mass = mass-mass1;
526    }
527    
528    fMinMass = mass-range;
529    fMaxMass = mass+range;
530    
531    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
532    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
533 }
534
535 //_______________________________________________________________________________
536
537 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
538    if(uplimit>lowlimit)
539    {
540       fMinMass = lowlimit;
541       fMaxMass = uplimit;
542    }
543    else{
544       printf("Error! Lower limit larger than upper limit!\n");
545       fMinMass = uplimit - uplimit*0.2;
546       fMaxMass = uplimit;
547    }
548 }
549
550 //_______________________________________________________________________________
551
552 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
553    if(nptbins>30) {
554       AliInfo("Maximum number of bins allowed is 30!");
555       return kFALSE;
556    }
557    if(!width) return kFALSE;
558    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
559    return kTRUE;
560 }
561
562 //_______________________________________________________________________________
563
564 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
565    if(!part || !jet){
566       printf("Particle or jet do not exist!\n");
567       return -99;
568    }
569    Double_t p[3],pj[3];
570    Bool_t okpp=part->PxPyPz(p);
571    Bool_t okpj=jet->PxPyPz(pj);
572    if(!okpp || !okpj){
573       printf("Problems getting momenta\n");
574       return -999;
575    }
576    Double_t pjet=jet->P();
577    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
578    return z;
579 }
580
581 //_______________________________________________________________________________
582
583 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
584    
585    // Statistics 
586    TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
587    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
588    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
589    hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
590    hstat->GetXaxis()->SetBinLabel(4,"N jets");
591    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
592    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
593    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
594    hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
595    hstat->SetNdivisions(1);
596    fmyOutput->Add(hstat);
597    
598    const Int_t nbinsmass=200;
599    const Int_t nbinsptjet=500;
600    const Int_t nbinsptD=100;
601    const Int_t nbinsz=100;
602    const Int_t nbinsphi=200;
603    const Int_t nbinseta=100;
604    
605    const Float_t ptjetlims[2]={0.,200.};
606    const Float_t ptDlims[2]={0.,50.};
607    const Float_t zlims[2]={0.,1.2};
608    const Float_t philims[2]={0.,6.3};
609    const Float_t etalims[2]={-1.5,1.5};
610    
611    if(fCandidateType==kDstartoKpipi) 
612    {
613       
614       TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
615       hDiffSideBand->SetStats(kTRUE);
616       hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
617       hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
618       hDiffSideBand->Sumw2();
619       fmyOutput->Add(hDiffSideBand); 
620       
621       
622       TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
623       hPtPion->SetStats(kTRUE);
624       hPtPion->GetXaxis()->SetTitle("GeV/c");
625       hPtPion->GetYaxis()->SetTitle("Entries");
626       hPtPion->Sumw2();
627       fmyOutput->Add(hPtPion);
628       
629    }
630    // D related histograms
631    TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
632    hNDPerEvNoJet->Sumw2();
633    fmyOutput->Add(hNDPerEvNoJet);
634    
635    TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
636    hptDPerEvNoJet->Sumw2();
637    fmyOutput->Add(hptDPerEvNoJet);
638    
639    const Int_t    nAxisD=4;
640    const Int_t    nbinsSparseD[nAxisD]={nbinsphi,nbinseta,nbinsptD,nbinsmass};
641    const Double_t minSparseD[nAxisD]  ={philims[0],etalims[0],ptDlims[0],fMinMass};
642    const Double_t maxSparseD[nAxisD]  ={philims[1],etalims[1],ptDlims[1],fMaxMass};
643    THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
644    hsDstandalone->Sumw2();
645    
646    fmyOutput->Add(hsDstandalone);
647    
648    // jet related fistograms
649    
650    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
651    hEjetTrks->Sumw2();
652    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
653    hPhiJetTrks->Sumw2();
654    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
655    hEtaJetTrks->Sumw2();
656    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
657    hPtJetTrks->Sumw2();
658    
659    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
660    hEjet->Sumw2();
661    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
662    hPhiJet->Sumw2();
663    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
664    hEtaJet->Sumw2();
665    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
666    hPtJet->Sumw2();
667    
668    TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
669    hPtJetWithD->Sumw2();
670    //for the MC this histogram is filled with the real background
671    TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
672    hPtJetWithDsb->Sumw2();
673    
674    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
675    hdeltaRJetTracks->Sumw2();
676    
677    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
678    hNtrArr->Sumw2();
679    
680    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
681    hNJetPerEv->Sumw2();
682    
683    TH1F *hNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
684    hNJetPerEvNoD->Sumw2();
685    
686    TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
687    hPtJetPerEvNoD->Sumw2();
688     
689    fmyOutput->Add(hEjetTrks);
690    fmyOutput->Add(hPhiJetTrks);
691    fmyOutput->Add(hEtaJetTrks);
692    fmyOutput->Add(hPtJetTrks);
693    fmyOutput->Add(hEjet);
694    fmyOutput->Add(hPhiJet);
695    fmyOutput->Add(hEtaJet);
696    fmyOutput->Add(hPtJet);
697    fmyOutput->Add(hPtJetWithD);
698    fmyOutput->Add(hPtJetWithDsb);
699    fmyOutput->Add(hdeltaRJetTracks);
700    fmyOutput->Add(hNtrArr);
701    fmyOutput->Add(hNJetPerEv);
702    fmyOutput->Add(hNJetPerEvNoD);
703    fmyOutput->Add(hPtJetPerEvNoD);
704    
705    TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
706    hDeltaRD->Sumw2();
707    fmyOutput->Add(hDeltaRD);
708    
709    //background (side bands for the Dstar and like sign for D0)
710    fJetRadius=GetJetContainer(0)->GetJetRadius();
711    TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
712    hInvMassptD->SetStats(kTRUE);
713    hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
714    hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
715    hInvMassptD->Sumw2();
716    
717    fmyOutput->Add(hInvMassptD);
718    
719    if(fUseMCInfo){
720       TH2F* hInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
721       hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
722       hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
723       hInvMassptDbg->Sumw2();
724       fmyOutput->Add(hInvMassptDbg);
725       
726    }
727    Double_t pi=TMath::Pi(), philims2[2];
728    philims2[0]=-pi*1./2.;
729    philims2[1]=pi*3./2.;
730    const Int_t nAxis=6;   
731    const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2};
732    const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5};
733    const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5};
734    THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass and IsSB", nAxis, nbinsSparse, minSparse, maxSparse);
735    hsDphiz->Sumw2();
736    
737    fmyOutput->Add(hsDphiz);
738
739    PostData(1, fmyOutput);
740    
741    return kTRUE; 
742 }
743
744 //_______________________________________________________________________________
745
746 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
747    
748    Double_t ptD=candidate->Pt();
749    Double_t ptjet=jet->Pt();
750    Double_t deltaR=DeltaR(candidate,jet);
751    Double_t deltaphi = jet->Phi()-candidate->Phi();
752    if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
753    if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
754    Double_t z=Z(candidate,jet);
755    
756    TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
757    hDeltaRD->Fill(deltaR); 
758    if(fUseReco){
759       if(fCandidateType==kD0toKpi) {
760          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
761          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
762          
763       }
764       
765       if(fCandidateType==kDstartoKpipi) {
766          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
767          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
768          
769       }
770    } else{ //generated level
771       //AliInfo("Non reco");
772       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
773    }
774 }
775
776 //_______________________________________________________________________________
777
778 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
779
780
781    Double_t masses[2]={0.,0.};
782    Int_t pdgdaughtersD0[2]={211,321};//pi,K 
783    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
784    
785    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
786    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
787    
788    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
789    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
790    Double_t point[6]={z,dPhi,ptj,ptD,masses[0],0};
791    
792    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
793    if(isselected==1 || isselected==3) {
794       
795       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
796       
797       FillMassHistograms(masses[0], ptD, deltaR);
798       hsDphiz->Fill(point,1.);
799    }
800    if(isselected>=2) {
801       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
802       
803       FillMassHistograms(masses[1], ptD, deltaR);
804       point[4]=masses[1];
805       hsDphiz->Fill(point,1.);
806    }
807    
808 }
809
810 //_______________________________________________________________________________
811
812 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
813   //dPhi and z not used at the moment,but will be (re)added
814
815    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
816    Double_t deltamass= dstar->DeltaInvMass(); 
817    //Double_t massD0= dstar->InvMassD0();
818    
819    
820    TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
821    hPtPion->Fill(softpi->Pt());
822    
823    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
824    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
825    Int_t isSB=IsDzeroSideBand(dstar);
826    
827    Double_t point[6]={z,dPhi,ptj,ptD,deltamass,isSB};
828
829    if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
830    
831    FillMassHistograms(deltamass, ptD, deltaR);
832    hsDphiz->Fill(point,1.);
833    
834 }
835
836 //_______________________________________________________________________________
837
838 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
839    
840    Double_t pdgmass=0;
841    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
842    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
843    Double_t point[6]={z,dPhi,ptjet,ptD,pdgmass,0};
844
845    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
846    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
847    point[4]=pdgmass;
848
849    if(deltaR<fJetRadius) {
850      hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
851      hsDphiz->Fill(point,1.);
852    }
853    
854 }
855
856 //_______________________________________________________________________________
857
858 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
859    
860    if(deltaR<fJetRadius) {
861       TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
862       hInvMassptD->Fill(mass,ptD);
863    }
864 }
865
866 //________________________________________________________________________________
867
868 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
869    Double_t deltaR=DeltaR(charm, jet);
870    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
871    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
872    if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
873    
874    return;
875    
876 }
877
878 //_______________________________________________________________________________
879
880 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
881    
882    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
883    // (expected detector resolution) on the left and right frm the D0 mass. Each band
884    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
885    TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
886    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
887    
888    Double_t deltaM=candDstar->DeltaInvMass(); 
889    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
890    Double_t ptD=candDstar->Pt();
891    Double_t ptjet=jet->Pt();
892    Double_t dPhi=jet->Phi()-candDstar->Phi();
893    Double_t deltaR=DeltaR(candDstar,jet);
894    if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
895    if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
896    
897    Int_t isSideBand=IsDzeroSideBand(candDstar);
898    //fill the background histogram with the side bands or when looking at MC with the real background
899    if(isSideBand==1){
900       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    
901       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
902       
903       if(deltaR<fJetRadius){  // evaluate in the near side      
904          //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
905          hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
906       }
907    }
908 }
909
910 //_______________________________________________________________________________
911
912 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
913    
914    //need mass, deltaR, pt jet, ptD
915    //two cases: D0 and Dstar
916    
917    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
918    if(!isselected) return;
919    
920    Double_t ptD=candbg->Pt();
921    Double_t ptjet=jet->Pt();
922    Double_t deltaR=DeltaR(candbg,jet);
923    
924    TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
925    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
926    
927    
928    if(fCandidateType==kDstartoKpipi){
929       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
930       Double_t deltaM=dstarbg->DeltaInvMass();
931       hInvMassptDbg->Fill(deltaM,ptD);
932       if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
933    }
934    
935    if(fCandidateType==kD0toKpi){
936       Double_t masses[2]={0.,0.};
937       Int_t pdgdaughtersD0[2]={211,321};//pi,K 
938       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
939       
940       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
941       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
942       
943       
944       if(isselected==1 || isselected==3) {
945          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
946          hInvMassptDbg->Fill(masses[0],ptD);
947       }
948       if(isselected>=2) {
949          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
950          hInvMassptDbg->Fill(masses[1],ptD);
951       }
952       
953       
954    }
955 }
956
957 //_______________________________________________________________________________
958
959 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
960    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
961    
962    if(!p1 || !p2) return -1;
963    Double_t phi1=p1->Phi(),eta1=p1->Eta();
964    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
965    
966    Double_t dPhi=phi1-phi2;
967    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
968    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
969    
970    Double_t dEta=eta1-eta2;
971    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
972    return deltaR;
973    
974 }
975
976 //_______________________________________________________________________________
977
978 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
979    
980    Double_t ptD=candDstar->Pt();
981    Int_t bin = fCuts->PtBin(ptD);
982    if (bin < 0){
983       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
984       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?      
985       return -1;  // out of bounds
986    }
987    
988    Double_t invM=candDstar->InvMassD0();
989    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
990
991    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
992    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
993    
994    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
995    else return 0;   
996
997 }