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