MCmodif1 EMCal phys sel off, EMCal trigger on
[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*)fSideBandArray->At(ib);
609                if(!charmbg) continue;
610                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
611                if (fIsDInJet) FlagFlavour(jet); //this are backgroud HF jets, but flagged as signal at the moment. Can use the bkg flavour flag in the future. This info is not stored now a part in the jet
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,jet);
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             if(fUseMCInfo){
1031                AliInfo("Creating a 7 axes container (MB background candidates)");
1032                const Int_t nAxis=9;   
1033                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1034                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1035                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1036                fNAxesBigSparse=nAxis;
1037                hsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1038                
1039             } else{
1040                AliInfo("Creating a 6 axes container");
1041                const Int_t nAxis=6;
1042                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1043                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1044                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1045                fNAxesBigSparse=nAxis;            
1046                
1047                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);
1048             }
1049          }
1050          if(!hsDphiz) AliFatal("No THnSparse created");
1051          hsDphiz->Sumw2();
1052          
1053          fOutput->Add(hsDphiz);
1054       }
1055    }
1056    PostData(1, fOutput);
1057    
1058    return kTRUE; 
1059 }
1060
1061 //_______________________________________________________________________________
1062
1063 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1064    
1065    Double_t ptD=candidate->Pt();
1066    Double_t deltaR=DeltaR(candidate,jet);
1067    Double_t phiD=candidate->Phi();
1068    Double_t deltaphi = jet->Phi()-phiD;
1069    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1070    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1071    Double_t z=Z(candidate,jet);
1072    /*
1073    if(z>1) {
1074       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1075       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1076       
1077       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1078       
1079       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1080       Double_t ptdiff=fPtJet-jet->Pt();
1081       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1082       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1083
1084       
1085    }
1086    */
1087    if(fIsDInJet)((TH1F*)fOutput->FindObject("hzDT"))->Fill(Z(candidate,jet,kTRUE));
1088    
1089    TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
1090    hDeltaRD->Fill(deltaR);
1091    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1092    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1093    if(fUseReco){
1094       if(fCandidateType==kD0toKpi) {
1095          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1096          
1097          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1098          
1099       }
1100       
1101       if(fCandidateType==kDstartoKpipi) {
1102          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1103          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1104          
1105       }
1106    } else{ //generated level
1107       //AliInfo("Non reco");
1108       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1109    }
1110 }
1111
1112 //_______________________________________________________________________________
1113
1114 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){
1115
1116
1117    Double_t masses[2]={0.,0.};
1118    Int_t pdgdaughtersD0[2]={211,321};//pi,K 
1119    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
1120    
1121    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1122    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1123    
1124    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1125    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1126    Double_t *point=0x0;
1127    if(fNAxesBigSparse==9){
1128       point=new Double_t[9];
1129       point[0]=z;
1130       point[1]=dPhi;
1131       point[2]=ptj;
1132       point[3]=ptD;
1133       point[4]=masses[0];
1134       point[5]=0;
1135       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1136       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1137       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1138    }
1139    if(fNAxesBigSparse==6){
1140       point=new Double_t[6];
1141       point[0]=z;
1142       point[1]=ptj;
1143       point[2]=ptD;
1144       point[3]=masses[0];
1145       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1146       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1147 }
1148   if(fNAxesBigSparse==7){
1149       point=new Double_t[7];
1150       point[0]=z;
1151       point[1]=ptj;
1152       point[2]=ptD;
1153       point[3]=masses[0];
1154       point[4]=0;
1155       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1156       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1157
1158    }
1159    
1160    
1161    Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());   
1162    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1163    if(isselected==1 || isselected==3) {
1164       
1165       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1166       
1167       FillMassHistograms(masses[0], ptD);
1168       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1169    }
1170    if(isselected>=2) {
1171       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1172       
1173       FillMassHistograms(masses[1], ptD);
1174       if(fNAxesBigSparse==9) point[4]=masses[1];
1175       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1176       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1177    }
1178    
1179 }
1180
1181 //_______________________________________________________________________________
1182
1183 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1184   //dPhi and z not used at the moment,but will be (re)added
1185
1186    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1187    Double_t deltamass= dstar->DeltaInvMass(); 
1188    //Double_t massD0= dstar->InvMassD0();
1189    
1190    
1191    TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
1192    hPtPion->Fill(softpi->Pt());
1193    
1194    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1195    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1196    Int_t isSB=0;//IsDzeroSideBand(dstar);
1197    
1198    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1199    Double_t *point=0x0;
1200    if(fNAxesBigSparse==9){
1201       point=new Double_t[9];
1202       point[0]=z;
1203       point[1]=dPhi;
1204       point[2]=ptj;
1205       point[3]=ptD;
1206       point[4]=deltamass;
1207       point[5]=static_cast<Double_t>(isSB);
1208       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1209       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1210       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1211    }
1212    if(fNAxesBigSparse==6){
1213       point=new Double_t[6];
1214       point[0]=z;
1215       point[1]=ptj;
1216       point[2]=ptD;
1217       point[3]=deltamass;
1218       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1219       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1220    }
1221    if(fNAxesBigSparse==7){
1222       point=new Double_t[7];
1223       point[0]=z;
1224       point[1]=ptj;
1225       point[2]=ptD;
1226       point[3]=deltamass;
1227       point[4]=0;
1228       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1229       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1230    }
1231
1232    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1233    
1234    FillMassHistograms(deltamass, ptD);
1235    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1236    
1237 }
1238
1239 //_______________________________________________________________________________
1240
1241 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1242    
1243    Double_t pdgmass=0;
1244    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1245    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1246    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1247    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1248    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1249    Double_t *point=0x0;
1250    if(fNAxesBigSparse==9){
1251       point=new Double_t[9];
1252       point[0]=z;
1253       point[1]=dPhi;
1254       point[2]=ptjet;
1255       point[3]=ptD;
1256       point[4]=pdgmass;
1257       point[5]=0;
1258       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1259       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1260       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1261    }
1262    if(fNAxesBigSparse==6){
1263       point=new Double_t[6];
1264       point[0]=z;
1265       point[1]=ptjet;
1266       point[2]=ptD;
1267       point[3]=pdgmass;
1268       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1269       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1270    }
1271       if(fNAxesBigSparse==7){
1272       point=new Double_t[6];
1273       point[0]=z;
1274       point[1]=ptjet;
1275       point[2]=ptD;
1276       point[3]=pdgmass;
1277       point[4]=1;
1278       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1279       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1280    }
1281
1282
1283    
1284    if(fNAxesBigSparse==9) point[4]=pdgmass;
1285    if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
1286    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1287    //if(fIsDInJet) {
1288    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1289    //}
1290    
1291 }
1292
1293 //_______________________________________________________________________________
1294
1295 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1296    
1297    if(fIsDInJet) {
1298       TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
1299       hInvMassptD->Fill(mass,ptD);
1300    }
1301 }
1302
1303 //________________________________________________________________________________
1304
1305 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1306    
1307    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1308    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1309    if (fIsDInJet) jet->AddFlavourTag(tag);
1310    
1311    return;
1312    
1313 }
1314
1315 //_______________________________________________________________________________
1316
1317 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1318    
1319    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
1320    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1321    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
1322    TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1323    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1324    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1325    
1326    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);  
1327    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1328    
1329    Double_t deltaM=candDstar->DeltaInvMass(); 
1330    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1331    Double_t z=Z(candDstar,jet);
1332    Double_t ptD=candDstar->Pt();
1333
1334    Double_t dPhi=jet->Phi()-candDstar->Phi();
1335
1336    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1337    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1338    
1339    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1340    //if no SB analysis we should not enter here, so no need of checking the number of axes
1341    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)};
1342    //fill the background histogram with the side bands or when looking at MC with the real background
1343    if(isSideBand==1){
1344       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    
1345       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1346       if(fSwitchOnSparses) hsDphiz->Fill(point,1.);
1347       //if(fIsDInJet){  // evaluate in the near side    
1348       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1349       //}
1350    }
1351 }
1352
1353 //_______________________________________________________________________________
1354
1355 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
1356    
1357    //need mass, deltaR, pt jet, ptD
1358    //two cases: D0 and Dstar
1359    
1360    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1361    if(!isselected) return;
1362    
1363    Double_t ptD=candbg->Pt();
1364    Double_t phiD=candbg->Phi();
1365    Double_t deltaphi = jet->Phi()-phiD;
1366    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1367    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1368    Double_t z=Z(candbg,jet);
1369
1370    Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1371    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1372
1373    TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1374    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1375
1376    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1377    Double_t *point=0x0;
1378    if(fNAxesBigSparse==9){
1379       point=new Double_t[9];
1380       point[0]=z;
1381       point[1]=deltaphi;
1382       point[2]=fPtJet;
1383       point[3]=ptD;
1384       point[4]=-999; //set below
1385       point[5]=1;
1386       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1387       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1388       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1389    }
1390
1391    if(fNAxesBigSparse==7){
1392       point=new Double_t[7];
1393       point[0]=z;
1394       point[1]=fPtJet;
1395       point[2]=ptD;
1396       point[3]=-999; //set below
1397       point[4]=1;
1398       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1399       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1400    }
1401
1402    if(fCandidateType==kDstartoKpipi){
1403       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1404       Double_t deltaM=dstarbg->DeltaInvMass();
1405       hInvMassptDbg->Fill(deltaM,ptD);
1406       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1407       if(fNAxesBigSparse==9) point[4]=deltaM;
1408       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
1409       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);      
1410    }
1411    
1412    if(fCandidateType==kD0toKpi){
1413       Double_t masses[2]={0.,0.};
1414       Int_t pdgdaughtersD0[2]={211,321};//pi,K 
1415       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
1416       
1417       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1418       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1419       
1420       
1421       if(isselected==1 || isselected==3) {
1422          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1423          hInvMassptDbg->Fill(masses[0],ptD);
1424          if(fNAxesBigSparse==9) point[4]=masses[0];
1425          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
1426          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1427      }
1428       if(isselected>=2) {
1429          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1430          hInvMassptDbg->Fill(masses[1],ptD);
1431          if(fNAxesBigSparse==9) point[4]=masses[1];
1432          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1433          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1434          
1435       }
1436       
1437       
1438    }
1439 }
1440
1441 //_______________________________________________________________________________
1442
1443 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1444    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1445    
1446    if(!p1 || !p2) return -1;
1447    Double_t phi1=p1->Phi(),eta1=p1->Eta();
1448    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1449    
1450    Double_t dPhi=phi1-phi2;
1451    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1452    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1453    
1454    Double_t dEta=eta1-eta2;
1455    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1456    return deltaR;
1457    
1458 }
1459
1460 //_______________________________________________________________________________
1461
1462 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1463    
1464    Double_t ptD=candDstar->Pt();
1465    Int_t bin = fCuts->PtBin(ptD);
1466    if (bin < 0){
1467       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1468       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?      
1469       return -1;  // out of bounds
1470    }
1471    
1472    Double_t invM=candDstar->InvMassD0();
1473    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1474
1475    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1476    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1477    
1478    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1479    else return 0;   
1480
1481 }
1482
1483 //_______________________________________________________________________________
1484
1485 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1486    //returns true/false and the array of particles not found among jet constituents
1487    
1488    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1489    TH1I* hIDddaugh   =(TH1I*)fOutput->FindObject("hIDddaugh");
1490    TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1491    TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1492    
1493    Int_t daughtersID[3];
1494    Int_t ndaugh=0;
1495    Int_t dmatchedID[3]={0,0,0};
1496    Int_t countmatches=0;
1497    if (fCandidateType==kDstartoKpipi){
1498       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1499       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1500       daughtersID[0]=dzero->GetProngID(0);
1501       daughtersID[1]=dzero->GetProngID(1);
1502       daughtersID[2]=dstar->GetBachelor()->GetID();
1503       ndaugh=3;
1504      
1505       dtrks=new AliAODTrack*[3];
1506       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1507       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1508       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1509   
1510       //check
1511       if(fillH){
1512          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  hControlDInJ->Fill(6);
1513          
1514          hIDddaugh->Fill(daughtersID[0]);
1515          hIDddaugh->Fill(daughtersID[1]);
1516          hIDddaugh->Fill(daughtersID[2]);
1517          
1518       }
1519       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1520    }
1521    
1522    if (fCandidateType==kD0toKpi){
1523       daughtersID[0]=charm->GetProngID(0);
1524       daughtersID[1]=charm->GetProngID(1);
1525       ndaugh=2;
1526       if(fillH){
1527          hIDddaugh->Fill(daughtersID[0]);
1528          hIDddaugh->Fill(daughtersID[1]);
1529       }
1530       dtrks=new AliAODTrack*[2];
1531       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1532       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1533
1534    }
1535    
1536    const Int_t cndaugh=ndaugh;
1537    daughOutOfJetID=new Int_t[cndaugh];
1538
1539    Int_t nchtrk=thejet->GetNumberOfTracks();
1540    for(Int_t itk=0;itk<nchtrk;itk++){
1541       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1542       Int_t idtkjet=tkjet->GetID();
1543       if(fillH) hIDjetTracks->Fill(idtkjet);
1544       for(Int_t id=0;id<ndaugh;id++){
1545          if(idtkjet==daughtersID[id]) {
1546             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1547             countmatches++;
1548             
1549          }
1550          
1551          if(countmatches==ndaugh) break;
1552       }
1553    }
1554    //reverse: include in the array the ID of daughters not matching
1555
1556    for(Int_t id=0;id<ndaugh;id++){
1557       if(dmatchedID[id]==0) {
1558          daughOutOfJetID[id]=daughtersID[id];
1559          if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1560       }
1561       else daughOutOfJetID[id]=0;
1562    }
1563    if(fillH){
1564       if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1565       if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1566       if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1567    }
1568    if(ndaugh!=countmatches){
1569       return kFALSE;
1570    }
1571    
1572    return kTRUE;
1573 }
1574
1575 //_______________________________________________________________________________
1576
1577 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1578    
1579    //check the conditions type:
1580    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet) 
1581    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1582    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1583     
1584    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1585    TH1F* hDRdaughOut=(TH1F*)fOutput->FindObject("hDRdaughOut");
1586    
1587    fPmissing[0]=0; 
1588    fPmissing[1]=0;
1589    fPmissing[2]=0;
1590    
1591    Bool_t testDeltaR=kFALSE;
1592    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1593    
1594    Int_t* daughOutOfJet;
1595    AliAODTrack** charmDaugh;
1596    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1597    if(testDaugh && testDeltaR) {
1598       //Z of candidates with daughters in the jet
1599       ((TH1F*)fOutput->FindObject("hzDinjet"))->Fill(Z(charm,thejet));
1600       
1601    }
1602    if(!testDaugh && testDeltaR && fTypeDInJet==2){
1603       
1604       Int_t ndaugh=3;
1605       if(fCandidateType==kD0toKpi) ndaugh=2;
1606       Int_t nOut=ndaugh;
1607       
1608       for(Int_t id=0;id<ndaugh;id++){
1609          if(daughOutOfJet[id]!=0){ //non-matched daughter
1610             //get track and its momentum
1611             nOut--;
1612             if(fillH) {
1613                hControlDInJ->Fill(3);
1614                if(id<2) hControlDInJ->Fill(4);
1615                if(id==2)hControlDInJ->Fill(5);
1616                hDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1617             }
1618             fPmissing[0]+=charmDaugh[id]->Px(); 
1619             fPmissing[1]+=charmDaugh[id]->Py();
1620             fPmissing[2]+=charmDaugh[id]->Pz();
1621          }
1622       
1623       }
1624       
1625       //now the D in within the jet
1626       testDaugh=kTRUE;
1627    
1628    }
1629    
1630    delete[] charmDaugh;
1631    
1632    Bool_t result=0;
1633    switch(fTypeDInJet){
1634    case 0:
1635       result=testDeltaR;
1636       break;
1637    case 1:
1638       result=testDeltaR && testDaugh;
1639       break;
1640    case 2:
1641       result=testDeltaR && testDaugh;
1642       break;
1643    default:
1644       AliInfo("Selection type not specified, use 1");
1645       result=testDeltaR && testDaugh;
1646       break;
1647    }
1648  return result;
1649 }
1650
1651 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1652    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1653    
1654    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1655    Bool_t binEMCal=kTRUE;
1656    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1657    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1658    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1659    return binEMCal;
1660
1661
1662 }