bd03629d4600a4b32aa211ede7212e386cba950d
[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 #include <TSystem.h>
34 #include <TObjectTable.h>
35
36 #include "AliAnalysisTaskFlavourJetCorrelations.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODHandler.h"
39 #include "AliAnalysisManager.h"
40 #include "AliLog.h"
41 #include "AliEmcalJet.h"
42 #include "AliJetContainer.h"
43 #include "AliAODRecoDecay.h"
44 #include "AliAODRecoCascadeHF.h"
45 #include "AliAODRecoDecayHF2Prong.h"
46 #include "AliESDtrack.h"
47 #include "AliAODMCParticle.h"
48 #include "AliPicoTrack.h"
49 #include "AliRDHFCutsD0toKpi.h"
50 #include "AliRDHFCutsDStartoKpipi.h"
51 #include "AliRhoParameter.h"
52
53 ClassImp(AliAnalysisTaskFlavourJetCorrelations)
54
55
56 //_______________________________________________________________________________
57
58 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
59 AliAnalysisTaskEmcalJet("",kTRUE),
60 fUseMCInfo(kTRUE), 
61 fUseReco(kTRUE),
62 fCandidateType(),
63 fPDGmother(),
64 fNProngs(),
65 fPDGdaughters(),
66 fBranchName(),
67 fCuts(0),
68 fMinMass(),
69 fMaxMass(),  
70 fJetArrName(0),
71 fCandArrName(0),
72 fLeadingJetOnly(kFALSE),
73 fJetRadius(0),
74 fCandidateArray(0),
75 fSideBandArray(0),
76 fJetOnlyMode(0),
77 fPmissing(),
78 fPtJet(0),
79 fIsDInJet(0),
80 fTypeDInJet(2),
81 fTrackArr(0),
82 fSwitchOnSB(0),
83 fSwitchOnPhiAxis(0),
84 fSwitchOnOutOfConeAxis(0),
85 fSwitchOnSparses(1),
86 fNAxesBigSparse(9)
87 {
88    //
89    // Default ctor
90    //
91    //SetMakeGeneralHistograms(kTRUE);
92    
93 }
94
95 //_______________________________________________________________________________
96
97 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, Bool_t jetOnly) :
98 AliAnalysisTaskEmcalJet(name,kTRUE),
99 fUseMCInfo(kTRUE),
100 fUseReco(kTRUE),  
101 fCandidateType(),
102 fPDGmother(),
103 fNProngs(),
104 fPDGdaughters(),
105 fBranchName(),
106 fCuts(0),
107 fMinMass(),
108 fMaxMass(),  
109 fJetArrName(0),
110 fCandArrName(0),
111 fLeadingJetOnly(kFALSE),
112 fJetRadius(0),
113 fCandidateArray(0),
114 fSideBandArray(0),
115 fJetOnlyMode(jetOnly),
116 fPmissing(),
117 fPtJet(0),
118 fIsDInJet(0),
119 fTypeDInJet(2),
120 fTrackArr(0),
121 fSwitchOnSB(0),
122 fSwitchOnPhiAxis(0),
123 fSwitchOnOutOfConeAxis(0),
124 fSwitchOnSparses(1),
125 fNAxesBigSparse(9)
126 {
127    //
128    // Constructor. Initialization of Inputs and Outputs
129    //
130    
131    Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
132    fCuts=cuts;
133    fCandidateType=candtype;
134    const Int_t nptbins=fCuts->GetNPtBins();
135    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};
136    
137    switch(fCandidateType){
138    case 0:
139       fPDGmother=421;
140       fNProngs=2;
141       fPDGdaughters[0]=211;//pi 
142       fPDGdaughters[1]=321;//K
143       fPDGdaughters[2]=0; //empty
144       fPDGdaughters[3]=0; //empty
145       fBranchName="D0toKpi";
146       fCandArrName="D0";
147       break;
148    case 1: 
149       fPDGmother=413;
150       fNProngs=3;
151       fPDGdaughters[1]=211;//pi soft
152       fPDGdaughters[0]=421;//D0
153       fPDGdaughters[2]=211;//pi fromD0
154       fPDGdaughters[3]=321; // K from D0
155       fBranchName="Dstar";
156       fCandArrName="Dstar";
157       
158       if(nptbins<=13){
159          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
160       } else {
161          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
162       }
163       break;
164    default:
165       printf("%d not accepted!!\n",fCandidateType);
166       break;
167    }
168    
169    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
170    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
171    
172    if(fJetOnlyMode){
173       DefineOutput(1,TList::Class()); //histos with jet info
174       DefineOutput(2,AliRDHFCuts::Class()); // my cuts
175    }
176    else{
177       DefineInput(1, TClonesArray::Class());
178       DefineInput(2, TClonesArray::Class());
179       
180       DefineOutput(1,TList::Class()); // histos
181       DefineOutput(2,AliRDHFCuts::Class()); // my cuts
182    }
183 }
184
185 //_______________________________________________________________________________
186
187 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
188    //
189    // destructor
190    //
191    
192    Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");  
193    
194    delete fCuts;
195    
196 }
197
198 //_______________________________________________________________________________
199
200 void AliAnalysisTaskFlavourJetCorrelations::Init(){
201    //
202    // Initialization
203    //
204    
205    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
206    
207    switch(fCandidateType){
208    case 0:
209       {
210          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
211          copyfCuts->SetName("AnalysisCutsDzero");
212          // Post the data
213          PostData(2,copyfCuts);
214       }
215       break;
216    case 1:
217       {
218          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
219          copyfCuts->SetName("AnalysisCutsDStar");
220          // Post the cuts
221          PostData(2,copyfCuts);
222       }
223       break;
224    default:
225       return;
226    }
227    
228    return;
229 }
230
231 //_______________________________________________________________________________
232
233 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { 
234    // output 
235    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
236    AliAnalysisTaskEmcal::UserCreateOutputObjects();
237    // define histograms
238    // the TList fOutput is already defined in  AliAnalysisTaskEmcal::UserCreateOutputObjects()
239    DefineHistoForAnalysis();
240    PostData(1,fOutput);
241    
242    return;
243 }
244
245 //_______________________________________________________________________________
246
247 Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
248 {
249    // user exec from AliAnalysisTaskEmcal is used
250     
251    // Load the event
252    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
253    
254    TClonesArray *arrayDStartoD0pi=0;
255    
256    if (!aodEvent && AODEvent() && IsStandardAOD()) {
257       
258       // In case there is an AOD handler writing a standard AOD, use the AOD 
259       // event in memory rather than the input (ESD) event.    
260       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
261       
262       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
263       // have to taken from the AOD event hold by the AliAODExtension
264       AliAODHandler* aodHandler = (AliAODHandler*) 
265       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
266       if(aodHandler->GetExtensions()) {
267          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
268          AliAODEvent *aodFromExt = ext->GetAOD();
269          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
270       }
271    } else {
272       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
273    }
274    
275    if (!arrayDStartoD0pi) {
276       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
277       //  return;
278    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
279    
280    TClonesArray* mcArray = 0x0;
281    if (fUseMCInfo) {
282       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
283       if (!mcArray) {
284          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
285          return kFALSE;
286       }
287    }
288    
289    //retrieve jets
290    fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
291    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
292    //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
293    fJetRadius=GetJetContainer(0)->GetJetRadius();
294    
295    if(!fTrackArr){
296       AliInfo(Form("Could not find the track array\n"));
297       return kFALSE;
298    }
299    
300     
301    fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
302    if (!fCandidateArray) return kFALSE;
303    if (fCandidateType==1 && fSwitchOnSB) {
304       fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
305       if (!fSideBandArray) return kFALSE;
306    }
307    //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
308    
309    //Histograms
310    TH1I* hstat=(TH1I*)fOutput->FindObject("hstat");
311    TH1F* hPtJetTrks=(TH1F*)fOutput->FindObject("hPtJetTrks");
312    TH1F* hPhiJetTrks=(TH1F*)fOutput->FindObject("hPhiJetTrks");
313    TH1F* hEtaJetTrks=(TH1F*)fOutput->FindObject("hEtaJetTrks");
314    TH1F* hEjetTrks=(TH1F*)fOutput->FindObject("hEjetTrks");
315    TH1F* hPtJet=(TH1F*)fOutput->FindObject("hPtJet");
316    TH1F* hPhiJet=(TH1F*)fOutput->FindObject("hPhiJet");
317    TH1F* hEtaJet=(TH1F*)fOutput->FindObject("hEtaJet");
318    TH1F* hEjet=(TH1F*)fOutput->FindObject("hEjet");
319    TH1F* hNtrArr=(TH1F*)fOutput->FindObject("hNtrArr");
320    TH1F* hNJetPerEv=(TH1F*)fOutput->FindObject("hNJetPerEv");
321    TH1F* hdeltaRJetTracks=(TH1F*)fOutput->FindObject("hdeltaRJetTracks");
322    TH1F* hNDPerEvNoJet=(TH1F*)fOutput->FindObject("hNDPerEvNoJet");
323    TH1F* hptDPerEvNoJet=(TH1F*)fOutput->FindObject("hptDPerEvNoJet");
324    TH1F* hNJetPerEvNoD=(TH1F*)fOutput->FindObject("hNJetPerEvNoD");
325    TH1F* hPtJetPerEvNoD=(TH1F*)fOutput->FindObject("hPtJetPerEvNoD");
326    THnSparseF* hnspDstandalone=(THnSparseF*)fOutput->FindObject("hsDstandalone");
327    THnSparseF* hsJet=(THnSparseF*)fOutput->FindObject("hsJet");
328    
329    TH1F* hztracksinjet=(TH1F*)fOutput->FindObject("hztracksinjet");
330    TH1F* hDiffPtTrPtJzNok=(TH1F*)fOutput->FindObject("hDiffPtTrPtJzNok");
331    TH1F* hDiffPtTrPtJzok=(TH1F*)fOutput->FindObject("hDiffPtTrPtJzok");
332    TH1F* hDiffPzTrPtJzok=(TH1F*)fOutput->FindObject("hDiffPzTrPtJzok");
333    TH1F* hDiffPzTrPtJzNok=(TH1F*)fOutput->FindObject("hDiffPzTrPtJzNok");
334    TH1I* hNtrkjzNok=(TH1I*)fOutput->FindObject("hNtrkjzNok");
335    TH1F* hztracksinjetT=(TH1F*)fOutput->FindObject("hztracksinjetT");
336    
337    
338    hstat->Fill(0);
339    
340    // fix for temporary bug in ESDfilter 
341    // the AODs with null vertex pointer didn't pass the PhysSel
342    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
343    
344    //Event selection
345    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
346    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
347    if(!iseventselected) return kFALSE;
348    
349    hstat->Fill(1);
350
351    //retrieve charm candidates selected
352    Int_t candidates = 0;
353    Int_t njets=GetJetContainer()->GetNJets();
354    
355    if(!fJetOnlyMode) {
356       candidates = fCandidateArray->GetEntriesFast();
357   
358    //trigger on jets  
359    if(njets == 0) {
360       hstat->Fill(6, candidates);
361       hNDPerEvNoJet->Fill(candidates);
362       for(Int_t iD=0;iD<candidates;iD++){
363          AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
364          if(!cand) continue;
365          hptDPerEvNoJet->Fill(cand->Pt());
366       
367       }
368       return kFALSE;
369       
370    }
371    
372    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
373    
374    for(Int_t ic = 0; ic < candidates; ic++) {
375       
376       // D* candidates
377       AliAODRecoDecayHF* charm=0x0;
378       AliAODRecoCascadeHF* dstar=0x0;
379       
380       
381       charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
382       if(!charm) continue;
383       
384       if(fCandidateType==kDstartoKpipi){ 
385          dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
386          if(!dstar) continue;
387       }
388       
389       hstat->Fill(2);
390       
391       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
392       
393       if(fCandidateType==kDstartoKpipi) {
394          if(fUseReco){
395             Double_t deltamass= dstar->DeltaInvMass();
396             candsparse[3]=deltamass;
397          } else candsparse[3] = 0.145; //for generated
398          if(fSwitchOnSparses) hnspDstandalone->Fill(candsparse);
399       }
400       if(fCandidateType==kD0toKpi){
401          if(fUseReco){
402             AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
403             Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
404             //not a further selection,just to get the value of isselected for the D0
405             Double_t masses[2];
406             Int_t pdgdaughtersD0[2]={211,321};//pi,K 
407             Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
408             
409             masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
410             masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
411             if(isselected==1 || isselected==3) {
412                candsparse[3]=masses[0];
413                if(fSwitchOnSparses) hnspDstandalone->Fill(candsparse);
414             }
415             if(isselected>=2){
416                candsparse[3]=masses[1];
417                if(fSwitchOnSparses) hnspDstandalone->Fill(candsparse);
418                
419             }
420          } else { //generated
421             Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
422             candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
423             if(fSwitchOnSparses) hnspDstandalone->Fill(candsparse);
424          }
425       }
426    }
427    }
428     
429     //Background Subtraction for jets
430     AliRhoParameter *rho = 0;
431     Double_t rhoval = 0;
432     TString sname("Rho");
433     if (!sname.IsNull()) {
434         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
435         if(rho) rhoval = rho->GetVal();
436     }
437
438    
439    // we start with jets
440    Double_t ejet   = 0;
441    Double_t phiJet = 0;
442    Double_t etaJet = 0;
443    Double_t ptjet = 0;
444    Double_t leadingJet =0;
445    Double_t pointJ[6];
446    
447    Int_t ntrarr=fTrackArr->GetEntriesFast();
448    hNtrArr->Fill(ntrarr);
449    
450    for(Int_t i=0;i<ntrarr;i++){
451       AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
452       //reusing histograms
453       hPtJetTrks->Fill(jtrack->Pt());
454       hPhiJetTrks->Fill(jtrack->Phi());
455       hEtaJetTrks->Fill(jtrack->Eta());
456       hEjetTrks->Fill(jtrack->E());
457    }
458    
459    
460    //option to use only the leading jet
461    if(fLeadingJetOnly){
462       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {    
463          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
464          ptjet   = jetL->Pt() - jetL->Area()*rhoval; //background subtraction
465          if(ptjet>leadingJet ) leadingJet = ptjet;
466          
467       }
468    }
469    
470    Int_t cntjet=0;
471    //loop over jets and consider only the leading jet in the event
472    for (Int_t iJets = 0; iJets<njets; iJets++) {
473       fPmissing[0]=0;
474       fPmissing[1]=0;
475       fPmissing[2]=0;
476       
477       //Printf("Jet N %d",iJets);
478       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
479       //Printf("Corr task Accept Jet");
480       if(!AcceptJet(jet)) {
481          hstat->Fill(5);
482          continue;
483       }
484       
485       //jets variables
486       ejet   = jet->E();
487       phiJet = jet->Phi();
488       etaJet = jet->Eta();
489       fPtJet = jet->Pt() - jet->Area()*rhoval; //background subtraction
490       Double_t origPtJet=fPtJet;
491       
492       pointJ[0] = phiJet;
493       pointJ[1] = etaJet;
494       pointJ[2] = ptjet - jet->Area()*rhoval; //background subtraction
495       pointJ[3] = ejet;
496       pointJ[4] = jet->GetNumberOfConstituents();
497       pointJ[5] = jet->Area();
498       
499       // choose the leading jet
500       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
501       //used for normalization
502       hstat->Fill(3);
503       cntjet++;
504       // fill energy, eta and phi of the jet
505       hEjet   ->Fill(ejet);
506       hPhiJet ->Fill(phiJet);
507       hEtaJet ->Fill(etaJet);
508       hPtJet  ->Fill(fPtJet);
509       if(fJetOnlyMode && fSwitchOnSparses) hsJet->Fill(pointJ,1);
510       //loop on jet particles
511       Int_t ntrjet=  jet->GetNumberOfTracks(); 
512       Double_t sumPtTracks=0,sumPzTracks=0;
513       Int_t zg1jtrk=0;
514       for(Int_t itrk=0;itrk<ntrjet;itrk++){
515          
516          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);     
517          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
518          Double_t ztrackj=Z(jetTrk,jet);
519          hztracksinjet->Fill(ztrackj);
520          sumPtTracks+=jetTrk->Pt(); 
521          sumPzTracks+=jetTrk->Pz(); 
522          if(ztrackj>1){
523             zg1jtrk++;
524          }
525          
526          Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
527          hztracksinjetT->Fill(ztrackjTr);
528          
529       }//end loop on jet tracks
530       
531       hNtrkjzNok->Fill(zg1jtrk);
532       Double_t diffpt=origPtJet-sumPtTracks;
533       Double_t diffpz=jet->Pz()-sumPzTracks;
534       if(zg1jtrk>0){
535          hDiffPtTrPtJzNok->Fill(diffpt);
536          hDiffPzTrPtJzNok->Fill(diffpz);
537       
538       }else{
539          hDiffPtTrPtJzok->Fill(diffpt);
540          hDiffPzTrPtJzok->Fill(diffpz);
541       }
542       
543       if(candidates==0){
544          hstat->Fill(7);
545          hPtJetPerEvNoD->Fill(fPtJet);
546       }
547       if(!fJetOnlyMode) {
548          //Printf("N candidates %d ", candidates);
549          for(Int_t ic = 0; ic < candidates; ic++) {
550             
551             // D* candidates
552             AliVParticle* charm=0x0;
553             charm=(AliVParticle*)fCandidateArray->At(ic);
554             if(!charm) continue;
555             AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
556             fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
557             if (fIsDInJet) FlagFlavour(jet);
558             if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
559             
560             //Note: the z component of the jet momentum comes from the eta-phi direction of the jet particles, it is not calculated from the z component of the tracks since, as default, the scheme used for jet reco is the pt-scheme which sums the scalar component, not the vectors. Addind the D daughter momentum component by componet as done here is not 100% correct, but the difference is small, for fairly collimated particles.
561
562             Double_t pjet[3];
563             jet->PxPyPz(pjet);
564              //background subtraction
565              pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
566              pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
567              pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
568             RecalculateMomentum(pjet,fPmissing);                            
569             fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
570             
571             
572             //debugging histograms
573             Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
574             for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
575             
576             ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
577             Double_t ptdiff=fPtJet-origPtJet;
578             ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
579             ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/origPtJet);
580             
581             FillHistogramsRecoJetCorr(charm, jet, aodEvent);
582             
583          }//end loop on candidates
584          
585          //retrieve side band background candidates for Dstar
586          Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
587          
588          for(Int_t ib=0;ib<nsbcand;ib++){
589             if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
590                AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
591                if(!sbcand) continue;
592                 
593                fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
594                Double_t pjet[3];
595                jet->PxPyPz(pjet);
596                 //background subtraction
597                 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
598                 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
599                 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
600                RecalculateMomentum(pjet,fPmissing);                                 
601                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
602                
603                SideBandBackground(sbcand,jet);
604                
605             }
606             if(fUseMCInfo){
607                AliAODRecoDecayHF* charmbg = 0x0;
608                charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
609                if(!charmbg) continue;
610                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
611                
612                Double_t pjet[3];
613                jet->PxPyPz(pjet);
614                 //background subtraction
615                 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
616                 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
617                 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
618                RecalculateMomentum(pjet,fPmissing);                                 
619                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
620                
621                MCBackground(charmbg);
622             }
623          }
624       }
625    } // end of jet loop
626    
627    hNJetPerEv->Fill(cntjet);
628    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
629    PostData(1,fOutput);
630    return kTRUE;
631 }
632
633 //_______________________________________________________________________________
634
635 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
636 {    
637    // The Terminate() function is the last function to be called during
638    // a query. It always runs on the client, it can be used to present
639    // the results graphically or save the results to file.
640    
641    Info("Terminate"," terminate");
642    AliAnalysisTaskSE::Terminate();
643    
644    fOutput = dynamic_cast<TList*> (GetOutputData(1));
645    if (!fOutput) {     
646       printf("ERROR: fOutput not available\n");
647       return;
648    }
649 }
650
651 //_______________________________________________________________________________
652
653 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
654    Float_t mass=0;
655    Int_t abspdg=TMath::Abs(pdg);
656    
657    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
658    // compute the Delta mass for the D*
659    if(fCandidateType==kDstartoKpipi){
660       Float_t mass1=0;
661       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
662       mass = mass-mass1;
663    }
664    
665    fMinMass = mass-range;
666    fMaxMass = mass+range;
667    
668    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
669    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
670 }
671
672 //_______________________________________________________________________________
673
674 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
675    if(uplimit>lowlimit)
676    {
677       fMinMass = lowlimit;
678       fMaxMass = uplimit;
679    }
680    else{
681       printf("Error! Lower limit larger than upper limit!\n");
682       fMinMass = uplimit - uplimit*0.2;
683       fMaxMass = uplimit;
684    }
685 }
686
687 //_______________________________________________________________________________
688
689 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
690    if(nptbins>30) {
691       AliInfo("Maximum number of bins allowed is 30!");
692       return kFALSE;
693    }
694    if(!width) return kFALSE;
695    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
696    return kTRUE;
697 }
698
699 //_______________________________________________________________________________
700
701 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
702    if(!part || !jet){
703       printf("Particle or jet do not exist!\n");
704       return -99;
705    }
706    Double_t p[3],pj[3];
707    Bool_t okpp=part->PxPyPz(p);
708    Bool_t okpj=jet->PxPyPz(pj);
709     
710     //Background Subtraction
711     AliRhoParameter *rho = 0;
712     Double_t rhoval = 0;
713     TString sname("Rho");
714     if (!sname.IsNull()) {
715         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
716         if(rho){
717             rhoval = rho->GetVal();
718             pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
719             pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
720             pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
721         }
722     }
723
724     
725     
726    if(!okpp || !okpj){
727       printf("Problems getting momenta\n");
728       return -999;
729    }
730    
731    RecalculateMomentum(pj, fPmissing);
732    if(transverse) return ZT(p,pj);
733    else
734    return Z(p,pj);
735 }
736
737 //_______________________________________________________________________________
738 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
739    
740    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
741    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
742    return z;
743 }
744
745
746 //_______________________________________________________________________________
747 Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
748    
749    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
750    Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
751    return z;
752 }
753
754 //_______________________________________________________________________________
755
756 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
757
758    pj[0]+=pmissing[0];
759    pj[1]+=pmissing[1];
760    pj[2]+=pmissing[2];
761
762 }
763
764 //_______________________________________________________________________________
765
766 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
767    
768    // Statistics 
769    TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
770    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
771    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
772    hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
773    hstat->GetXaxis()->SetBinLabel(4,"N jets");
774    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
775    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
776    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
777    hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
778    hstat->SetNdivisions(1);
779    fOutput->Add(hstat);
780    
781    const Int_t nbinsmass=200;
782    const Int_t nbinsptjet=500;
783    const Int_t nbinsptD=100;
784    const Int_t nbinsz=100;
785    const Int_t nbinsphi=200;
786    const Int_t nbinseta=100;
787    
788    //binning for THnSparse
789    const Int_t nbinsSpsmass=50;
790    const Int_t nbinsSpsptjet=100;
791    const Int_t nbinsSpsptD=50;
792    const Int_t nbinsSpsz=100;
793    const Int_t nbinsSpsphi=100;
794    const Int_t nbinsSpseta=60;
795    const Int_t nbinsSpsContrib=100;
796    const Int_t nbinsSpsA=100;
797     
798    const Float_t ptjetlims[2]={0.,200.};
799    const Float_t ptDlims[2]={0.,50.};
800    const Float_t zlims[2]={0.,1.2};
801    const Float_t philims[2]={0.,6.3};
802    const Float_t etalims[2]={-1.5,1.5};
803    const Int_t   nContriblims[2]={0,100};
804    const Float_t arealims[2]={0.,2};
805    
806    // jet related fistograms
807    
808    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
809    hEjetTrks->Sumw2();
810    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
811    hPhiJetTrks->Sumw2();
812    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
813    hEtaJetTrks->Sumw2();
814    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
815    hPtJetTrks->Sumw2();
816    
817    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
818    hEjet->Sumw2();
819    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
820    hPhiJet->Sumw2();
821    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
822    hEtaJet->Sumw2();
823    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
824    hPtJet->Sumw2();
825    
826    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
827    hdeltaRJetTracks->Sumw2();
828    
829    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
830    hNtrArr->Sumw2();
831    
832    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
833    hNJetPerEv->Sumw2();
834    
835    
836    fOutput->Add(hEjetTrks);
837    fOutput->Add(hPhiJetTrks);
838    fOutput->Add(hEtaJetTrks);
839    fOutput->Add(hPtJetTrks);
840    fOutput->Add(hEjet);
841    fOutput->Add(hPhiJet);
842    fOutput->Add(hEtaJet);
843    fOutput->Add(hPtJet);
844    fOutput->Add(hdeltaRJetTracks);
845    fOutput->Add(hNtrArr);
846    fOutput->Add(hNJetPerEv);
847    
848    if(fJetOnlyMode && fSwitchOnSparses){
849       //thnsparse for jets
850       const Int_t nAxis=6;   
851       const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
852       const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
853       const Double_t 
854         maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
855       THnSparseF *hsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
856       hsJet->Sumw2();
857       
858       fOutput->Add(hsJet);
859    
860    }
861
862    if(!fJetOnlyMode){
863       
864       //debugging histograms
865       TH1I* hControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
866       hControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
867       hControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
868       hControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
869       hControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
870       hControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
871       hControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
872       hControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
873       hControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
874       
875       hControlDInJ->SetNdivisions(1);
876       hControlDInJ->GetXaxis()->SetLabelSize(0.05);
877       fOutput->Add(hControlDInJ);
878       
879       TH1F *hmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
880       fOutput->Add(hmissingp);
881       
882       for(Int_t k=0;k<3;k++) {
883          TH1F *hMissPi=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
884          fOutput->Add(hMissPi);
885       }
886       TH1F *hDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
887       
888       fOutput->Add(hDeltaPtJet);
889       TH1F *hRelDeltaPtJet=new TH1F("hRelDeltaPtJet", "Difference between the jet pt before and after correction/ original pt;(p_{T}^{jet,recalc}-p_{T}^{jet,orig})/p_{T}^{jet,orig}", 200, 0.,1.);
890       fOutput->Add(hRelDeltaPtJet);
891       
892       TH1F* hzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
893       fOutput->Add(hzDinjet);
894       //frag func of tracks in the jet
895       TH1F* hztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
896       fOutput->Add(hztracksinjet);
897       
898       //check jet and tracks
899       TH1F* hDiffPtTrPtJzok=new TH1F("hDiffPtTrPtJzok","Sum p_{T}^{trks} - p_{T}^{jet}, for z<1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
900       fOutput->Add(hDiffPtTrPtJzok);
901       TH1F* hDiffPtTrPtJzNok=new TH1F("hDiffPtTrPtJzNok","Sum p_{T}^{trks} - p_{T}^{jet}, for z>1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
902       fOutput->Add(hDiffPtTrPtJzNok);
903       TH1F* hDiffPzTrPtJzok=new TH1F("hDiffPzTrPtJzok","Sum p_{z}^{trks} - p_{z}^{jet}, for z<1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
904       fOutput->Add(hDiffPzTrPtJzok);
905       TH1F* hDiffPzTrPtJzNok=new TH1F("hDiffPzTrPtJzNok","Sum p_{z}^{trks} - p_{z}^{jet}, for z>1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
906       fOutput->Add(hDiffPzTrPtJzNok);
907       TH1I* hNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
908       fOutput->Add(hNtrkjzNok);
909       
910       //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
911       TH1F* hzDT=new TH1F("hzDT", "Z of D in jet in transverse components;(p_{T}^{D} dot p_{T}^{jet})/p_{T}^{jet}^{2} ",nbinsz,zlims[0],zlims[1]);
912       fOutput->Add(hzDT);
913       TH1F* hztracksinjetT=new TH1F("hztracksinjetT", "Z of jet tracks in transverse components;(p_{T}^{trks} dot p_{T}^{jet})/p_{T}^{jet}^{2}",nbinsz,zlims[0],zlims[1]);
914       fOutput->Add(hztracksinjetT);
915       
916       TH1I* hIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
917       fOutput->Add(hIDddaugh);
918       TH1I* hIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
919       fOutput->Add(hIDddaughOut);
920       TH1I* hIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
921       fOutput->Add(hIDjetTracks);
922       
923       TH1F* hDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
924       fOutput->Add(hDRdaughOut);
925       
926       
927       if(fCandidateType==kDstartoKpipi) 
928       {
929          
930          TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
931          hDiffSideBand->SetStats(kTRUE);
932          hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
933          hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
934          hDiffSideBand->Sumw2();
935          fOutput->Add(hDiffSideBand); 
936          
937          
938          TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
939          hPtPion->SetStats(kTRUE);
940          hPtPion->GetXaxis()->SetTitle("GeV/c");
941          hPtPion->GetYaxis()->SetTitle("Entries");
942          hPtPion->Sumw2();
943          fOutput->Add(hPtPion);
944          
945       }
946       // D related histograms
947       TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
948       hNDPerEvNoJet->Sumw2();
949       fOutput->Add(hNDPerEvNoJet);
950       
951       TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
952       hptDPerEvNoJet->Sumw2();
953       fOutput->Add(hptDPerEvNoJet);
954       
955       if(fSwitchOnSparses){
956          const Int_t    nAxisD=4;
957          const Int_t    nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
958          const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
959          const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
960          THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
961          hsDstandalone->Sumw2();
962          
963          fOutput->Add(hsDstandalone);
964       }
965       
966       /*
967       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]);
968       hPtJetWithD->Sumw2();
969       //for the MC this histogram is filled with the real background
970       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]);
971       hPtJetWithDsb->Sumw2();
972       
973       fOutput->Add(hPtJetWithD);
974       fOutput->Add(hPtJetWithDsb);
975
976       */
977       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);
978       hNJetPerEvNoD->Sumw2();
979       
980       TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
981       hPtJetPerEvNoD->Sumw2();
982       
983       fOutput->Add(hNJetPerEvNoD);
984       fOutput->Add(hPtJetPerEvNoD);
985       
986       TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
987       hDeltaRD->Sumw2();
988       fOutput->Add(hDeltaRD);
989       
990       //background (side bands for the Dstar and like sign for D0)
991       fJetRadius=GetJetContainer(0)->GetJetRadius();
992       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]);
993       hInvMassptD->SetStats(kTRUE);
994       hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
995       hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
996       hInvMassptD->Sumw2();
997       
998       fOutput->Add(hInvMassptD);
999       
1000       if(fUseMCInfo){
1001          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]);
1002          hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1003          hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1004          hInvMassptDbg->Sumw2();
1005          fOutput->Add(hInvMassptDbg);
1006          
1007       }
1008       
1009       if(fSwitchOnSparses){
1010          Double_t pi=TMath::Pi(), philims2[2];
1011          philims2[0]=-pi*1./2.;
1012          philims2[1]=pi*3./2.;
1013          THnSparseF *hsDphiz=0x0; //definition below according to the switches
1014          
1015          if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1016             AliInfo("Creating a 9 axes container: This might cause large memory usage");
1017             const Int_t nAxis=9;   
1018             const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1019             const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1020             const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1021             fNAxesBigSparse=nAxis;
1022             
1023             hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. SB? D within the jet cone?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1024          }
1025          
1026          if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1027             fSwitchOnPhiAxis=0;
1028             fSwitchOnOutOfConeAxis=0;
1029             fSwitchOnSB=0;
1030             AliInfo("Creating a 6 axes container");
1031             const Int_t nAxis=6;
1032             const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1033             const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1034             const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1035             fNAxesBigSparse=nAxis;               
1036             
1037             hsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass., D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1038          }
1039          if(!hsDphiz) AliFatal("No THnSparse created");
1040          hsDphiz->Sumw2();
1041          
1042          fOutput->Add(hsDphiz);
1043       }
1044    }
1045    PostData(1, fOutput);
1046    
1047    return kTRUE; 
1048 }
1049
1050 //_______________________________________________________________________________
1051
1052 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1053    
1054    Double_t ptD=candidate->Pt();
1055    Double_t deltaR=DeltaR(candidate,jet);
1056    Double_t phiD=candidate->Phi();
1057    Double_t deltaphi = jet->Phi()-phiD;
1058    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1059    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1060    Double_t z=Z(candidate,jet);
1061    /*
1062    if(z>1) {
1063       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1064       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1065       
1066       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1067       
1068       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1069       Double_t ptdiff=fPtJet-jet->Pt();
1070       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1071       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1072
1073       
1074    }
1075    */
1076    if(fIsDInJet)((TH1F*)fOutput->FindObject("hzDT"))->Fill(Z(candidate,jet,kTRUE));
1077    
1078    TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
1079    hDeltaRD->Fill(deltaR);
1080    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1081    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1082    if(fUseReco){
1083       if(fCandidateType==kD0toKpi) {
1084          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1085          
1086          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1087          
1088       }
1089       
1090       if(fCandidateType==kDstartoKpipi) {
1091          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1092          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1093          
1094       }
1095    } else{ //generated level
1096       //AliInfo("Non reco");
1097       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1098    }
1099 }
1100
1101 //_______________________________________________________________________________
1102
1103 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent){
1104
1105
1106    Double_t masses[2]={0.,0.};
1107    Int_t pdgdaughtersD0[2]={211,321};//pi,K 
1108    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
1109    
1110    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1111    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1112    
1113    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1114    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1115    Double_t *point=0x0;
1116    if(fNAxesBigSparse==9){
1117       point=new Double_t[9];
1118       point[0]=z;
1119       point[1]=dPhi;
1120       point[2]=ptj;
1121       point[3]=ptD;
1122       point[4]=masses[0];
1123       point[5]=0;
1124       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1125       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1126       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1127    }
1128    if(fNAxesBigSparse==6){
1129       point=new Double_t[6];
1130       point[0]=z;
1131       point[1]=ptj;
1132       point[2]=ptD;
1133       point[3]=masses[0];
1134       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1135       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1136    }
1137    
1138    
1139    Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());   
1140    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1141    if(isselected==1 || isselected==3) {
1142       
1143       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1144       
1145       FillMassHistograms(masses[0], ptD);
1146       if(fSwitchOnOutOfConeAxis || fIsDInJet) hsDphiz->Fill(point,1.);
1147    }
1148    if(isselected>=2) {
1149       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1150       
1151       FillMassHistograms(masses[1], ptD);
1152       if(fNAxesBigSparse==9) point[4]=masses[1];
1153       if(fNAxesBigSparse==6) point[3]=masses[1];
1154       if(fSwitchOnOutOfConeAxis || fIsDInJet) hsDphiz->Fill(point,1.);
1155    }
1156    
1157 }
1158
1159 //_______________________________________________________________________________
1160
1161 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1162   //dPhi and z not used at the moment,but will be (re)added
1163
1164    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1165    Double_t deltamass= dstar->DeltaInvMass(); 
1166    //Double_t massD0= dstar->InvMassD0();
1167    
1168    
1169    TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
1170    hPtPion->Fill(softpi->Pt());
1171    
1172    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1173    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1174    Int_t isSB=0;//IsDzeroSideBand(dstar);
1175    
1176    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1177    Double_t *point=0x0;
1178    if(fNAxesBigSparse==9){
1179       point=new Double_t[9];
1180       point[0]=z;
1181       point[1]=dPhi;
1182       point[2]=ptj;
1183       point[3]=ptD;
1184       point[4]=deltamass;
1185       point[5]=static_cast<Double_t>(isSB);
1186       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1187       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1188       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1189    }
1190    if(fNAxesBigSparse==6){
1191       point=new Double_t[6];
1192       point[0]=z;
1193       point[1]=ptj;
1194       point[2]=ptD;
1195       point[3]=deltamass;
1196       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1197       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1198    }
1199
1200    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1201    
1202    FillMassHistograms(deltamass, ptD);
1203    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1204    
1205 }
1206
1207 //_______________________________________________________________________________
1208
1209 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1210    
1211    Double_t pdgmass=0;
1212    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1213    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1214    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1215    Double_t *point=0x0;
1216    if(fNAxesBigSparse==9){
1217       point=new Double_t[9];
1218       point[0]=z;
1219       point[1]=dPhi;
1220       point[2]=ptjet;
1221       point[3]=ptD;
1222       point[4]=pdgmass;
1223       point[5]=0;
1224       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1225       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1226       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1227    }
1228    if(fNAxesBigSparse==6){
1229       point=new Double_t[6];
1230       point[0]=z;
1231       point[1]=ptjet;
1232       point[2]=ptD;
1233       point[3]=pdgmass;
1234       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1235       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1236    }
1237
1238    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1239    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1240    if(fNAxesBigSparse==9) point[4]=pdgmass;
1241    if(fNAxesBigSparse==6) point[3]=pdgmass;
1242    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1243    //if(fIsDInJet) {
1244    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1245    //}
1246    
1247 }
1248
1249 //_______________________________________________________________________________
1250
1251 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1252    
1253    if(fIsDInJet) {
1254       TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
1255       hInvMassptD->Fill(mass,ptD);
1256    }
1257 }
1258
1259 //________________________________________________________________________________
1260
1261 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1262    
1263    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1264    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1265    if (fIsDInJet) jet->AddFlavourTag(tag);
1266    
1267    return;
1268    
1269 }
1270
1271 //_______________________________________________________________________________
1272
1273 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1274    
1275    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
1276    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1277    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
1278    TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1279    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1280    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1281    
1282    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);  
1283    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1284    
1285    Double_t deltaM=candDstar->DeltaInvMass(); 
1286    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1287    Double_t z=Z(candDstar,jet);
1288    Double_t ptD=candDstar->Pt();
1289
1290    Double_t dPhi=jet->Phi()-candDstar->Phi();
1291
1292    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1293    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1294    
1295    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1296    //if no SB analysis we should not enter here, so no need of checking the number of axes
1297    Double_t point[9]={z,dPhi,fPtJet,ptD,deltaM,static_cast<Double_t>(isSideBand), static_cast<Double_t>(fIsDInJet ? 1 : 0),static_cast<Double_t>(bDInEMCalAcc? 1 : 0),static_cast<Double_t>(bJetInEMCalAcc? 1 : 0)};
1298    //fill the background histogram with the side bands or when looking at MC with the real background
1299    if(isSideBand==1){
1300       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    
1301       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1302       if(fSwitchOnSparses) hsDphiz->Fill(point,1.);
1303       //if(fIsDInJet){  // evaluate in the near side    
1304       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1305       //}
1306    }
1307 }
1308
1309 //_______________________________________________________________________________
1310
1311 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg){
1312    
1313    //need mass, deltaR, pt jet, ptD
1314    //two cases: D0 and Dstar
1315    
1316    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1317    if(!isselected) return;
1318    
1319    Double_t ptD=candbg->Pt();
1320    
1321    TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1322    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1323
1324
1325    if(fCandidateType==kDstartoKpipi){
1326       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1327       Double_t deltaM=dstarbg->DeltaInvMass();
1328       hInvMassptDbg->Fill(deltaM,ptD);
1329       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1330    }
1331    
1332    if(fCandidateType==kD0toKpi){
1333       Double_t masses[2]={0.,0.};
1334       Int_t pdgdaughtersD0[2]={211,321};//pi,K 
1335       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
1336       
1337       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1338       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1339       
1340       
1341       if(isselected==1 || isselected==3) {
1342          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1343          hInvMassptDbg->Fill(masses[0],ptD);
1344       }
1345       if(isselected>=2) {
1346          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1347          hInvMassptDbg->Fill(masses[1],ptD);
1348       }
1349       
1350       
1351    }
1352 }
1353
1354 //_______________________________________________________________________________
1355
1356 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1357    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1358    
1359    if(!p1 || !p2) return -1;
1360    Double_t phi1=p1->Phi(),eta1=p1->Eta();
1361    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1362    
1363    Double_t dPhi=phi1-phi2;
1364    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1365    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1366    
1367    Double_t dEta=eta1-eta2;
1368    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1369    return deltaR;
1370    
1371 }
1372
1373 //_______________________________________________________________________________
1374
1375 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1376    
1377    Double_t ptD=candDstar->Pt();
1378    Int_t bin = fCuts->PtBin(ptD);
1379    if (bin < 0){
1380       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1381       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?      
1382       return -1;  // out of bounds
1383    }
1384    
1385    Double_t invM=candDstar->InvMassD0();
1386    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1387
1388    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1389    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1390    
1391    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1392    else return 0;   
1393
1394 }
1395
1396 //_______________________________________________________________________________
1397
1398 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1399    //returns true/false and the array of particles not found among jet constituents
1400    
1401    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1402    TH1I* hIDddaugh   =(TH1I*)fOutput->FindObject("hIDddaugh");
1403    TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1404    TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1405    
1406    Int_t daughtersID[3];
1407    Int_t ndaugh=0;
1408    Int_t dmatchedID[3]={0,0,0};
1409    Int_t countmatches=0;
1410    if (fCandidateType==kDstartoKpipi){
1411       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1412       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1413       daughtersID[0]=dzero->GetProngID(0);
1414       daughtersID[1]=dzero->GetProngID(1);
1415       daughtersID[2]=dstar->GetBachelor()->GetID();
1416       ndaugh=3;
1417      
1418       dtrks=new AliAODTrack*[3];
1419       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1420       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1421       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1422   
1423       //check
1424       if(fillH){
1425          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  hControlDInJ->Fill(6);
1426          
1427          hIDddaugh->Fill(daughtersID[0]);
1428          hIDddaugh->Fill(daughtersID[1]);
1429          hIDddaugh->Fill(daughtersID[2]);
1430          
1431       }
1432       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1433    }
1434    
1435    if (fCandidateType==kD0toKpi){
1436       daughtersID[0]=charm->GetProngID(0);
1437       daughtersID[1]=charm->GetProngID(1);
1438       ndaugh=2;
1439       if(fillH){
1440          hIDddaugh->Fill(daughtersID[0]);
1441          hIDddaugh->Fill(daughtersID[1]);
1442       }
1443       dtrks=new AliAODTrack*[2];
1444       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1445       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1446
1447    }
1448    
1449    const Int_t cndaugh=ndaugh;
1450    daughOutOfJetID=new Int_t[cndaugh];
1451
1452    Int_t nchtrk=thejet->GetNumberOfTracks();
1453    for(Int_t itk=0;itk<nchtrk;itk++){
1454       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1455       Int_t idtkjet=tkjet->GetID();
1456       if(fillH) hIDjetTracks->Fill(idtkjet);
1457       for(Int_t id=0;id<ndaugh;id++){
1458          if(idtkjet==daughtersID[id]) {
1459             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1460             countmatches++;
1461             
1462          }
1463          
1464          if(countmatches==ndaugh) break;
1465       }
1466    }
1467    //reverse: include in the array the ID of daughters not matching
1468
1469    for(Int_t id=0;id<ndaugh;id++){
1470       if(dmatchedID[id]==0) {
1471          daughOutOfJetID[id]=daughtersID[id];
1472          if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1473       }
1474       else daughOutOfJetID[id]=0;
1475    }
1476    if(fillH){
1477       if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1478       if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1479       if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1480    }
1481    if(ndaugh!=countmatches){
1482       return kFALSE;
1483    }
1484    
1485    return kTRUE;
1486 }
1487
1488 //_______________________________________________________________________________
1489
1490 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1491    
1492    //check the conditions type:
1493    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet) 
1494    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1495    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1496     
1497    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1498    TH1F* hDRdaughOut=(TH1F*)fOutput->FindObject("hDRdaughOut");
1499    
1500    fPmissing[0]=0; 
1501    fPmissing[1]=0;
1502    fPmissing[2]=0;
1503    
1504    Bool_t testDeltaR=kFALSE;
1505    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1506    
1507    Int_t* daughOutOfJet;
1508    AliAODTrack** charmDaugh;
1509    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1510    if(testDaugh && testDeltaR) {
1511       //Z of candidates with daughters in the jet
1512       ((TH1F*)fOutput->FindObject("hzDinjet"))->Fill(Z(charm,thejet));
1513       
1514    }
1515    if(!testDaugh && testDeltaR && fTypeDInJet==2){
1516       
1517       Int_t ndaugh=3;
1518       if(fCandidateType==kD0toKpi) ndaugh=2;
1519       Int_t nOut=ndaugh;
1520       
1521       for(Int_t id=0;id<ndaugh;id++){
1522          if(daughOutOfJet[id]!=0){ //non-matched daughter
1523             //get track and its momentum
1524             nOut--;
1525             if(fillH) {
1526                hControlDInJ->Fill(3);
1527                if(id<2) hControlDInJ->Fill(4);
1528                if(id==2)hControlDInJ->Fill(5);
1529                hDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1530             }
1531             fPmissing[0]+=charmDaugh[id]->Px(); 
1532             fPmissing[1]+=charmDaugh[id]->Py();
1533             fPmissing[2]+=charmDaugh[id]->Pz();
1534          }
1535       
1536       }
1537       
1538       //now the D in within the jet
1539       testDaugh=kTRUE;
1540    
1541    }
1542    
1543    delete[] charmDaugh;
1544    
1545    Bool_t result=0;
1546    switch(fTypeDInJet){
1547    case 0:
1548       result=testDeltaR;
1549       break;
1550    case 1:
1551       result=testDeltaR && testDaugh;
1552       break;
1553    case 2:
1554       result=testDeltaR && testDaugh;
1555       break;
1556    default:
1557       AliInfo("Selection type not specified, use 1");
1558       result=testDeltaR && testDaugh;
1559       break;
1560    }
1561  return result;
1562 }
1563
1564 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1565    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1566    
1567    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1568    Bool_t binEMCal=kTRUE;
1569    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1570    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1571    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1572    return binEMCal;
1573
1574
1575 }