Add ExchangeContainer and Fix bug Clear TClonesArray
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.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 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A.Grelli (Utrecht University) a.grelli@uu.nl
23 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
24 //-----------------------------------------------------------------------
25
26 #include <TDatabasePDG.h>
27 #include <TParticle.h>
28 #include <TVector3.h>
29 #include "TROOT.h"
30 #include <TH3F.h>
31
32 #include "AliRDHFCutsDStartoKpipi.h"
33 #include "AliRDHFCutsD0toKpi.h"
34 #include "AliAODMCHeader.h"
35 #include "AliAODHandler.h"
36 #include "AliAnalysisManager.h"
37 #include "AliLog.h"
38 #include "AliAODVertex.h"
39 #include "AliAODRecoDecay.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAODRecoDecayHF2Prong.h"
42 #include "AliESDtrack.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
45
46 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
47
48 //_______________________________________________________________________________
49
50 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
51 AliAnalysisTaskSE(),
52 fUseMCInfo(kFALSE),
53 fUseReco(kTRUE),
54 fCandidateType(0),
55 fCandidateName(""),
56 fPDGmother(0),
57 fNProngs(0),
58 fBranchName(""),
59 fOutput(0),
60 fCuts(0),
61 fMinMass(0.),
62 fMaxMass(0.),
63 fCandidateArray(0),
64 fSideBandArray(0)
65
66 {
67    //
68    // Default constructor
69    //
70    
71    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
72    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
73 }
74
75 //_______________________________________________________________________________
76
77 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
78 AliAnalysisTaskSE(name),
79 fUseMCInfo(kFALSE),
80 fUseReco(kTRUE),
81 fCandidateType(candtype),
82 fCandidateName(""),
83 fPDGmother(0),
84 fNProngs(0),
85 fBranchName(""),
86 fOutput(0),
87 fCuts(cuts),
88 fMinMass(0.),
89 fMaxMass(0.),
90 fCandidateArray(0),
91 fSideBandArray(0)
92 {
93    //
94    // Constructor. Initialization of Inputs and Outputs
95    //
96    
97    Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
98    
99    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
100    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
101    
102    const Int_t nptbins = fCuts->GetNPtBins();
103    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 };
104    
105    switch (fCandidateType) {
106    case 0 :
107       fCandidateName = "D0";
108       fPDGmother = 421;
109       fNProngs = 2;
110       fPDGdaughters[0] = 211;  // pi 
111       fPDGdaughters[1] = 321;  // K
112       fPDGdaughters[2] = 0;    // empty
113       fPDGdaughters[3] = 0;    // empty
114       fBranchName = "D0toKpi";
115       break;
116    case 1 :
117       fCandidateName = "Dstar";
118       fPDGmother = 413;
119       fNProngs = 3;
120       fPDGdaughters[1] = 211; // pi soft
121       fPDGdaughters[0] = 421; // D0
122       fPDGdaughters[2] = 211; // pi fromD0
123       fPDGdaughters[3] = 321; // K from D0
124       fBranchName = "Dstar";
125       
126       if (nptbins<=13) {
127          for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
128       } else {
129          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
130       }
131       break;
132    default :
133       printf("%d not accepted!!\n",fCandidateType);
134       break;
135    }
136    
137    if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
138    if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
139    
140    DefineOutput(1, TList::Class());       // histos
141    DefineOutput(2, AliRDHFCuts::Class()); // my cuts
142    DefineOutput(3, TClonesArray::Class()); //array of candidates
143    DefineOutput(4, TClonesArray::Class()); //array of SB candidates
144 }
145
146 //_______________________________________________________________________________
147
148 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
149 {
150    //
151    // destructor
152    //
153    
154    Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  
155    
156    if (fOutput) { delete fOutput; fOutput = 0; }
157    if (fCuts)   { delete fCuts;   fCuts   = 0; }
158    if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
159    delete fSideBandArray;
160    
161 }
162
163 //_______________________________________________________________________________
164
165 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
166 {
167    //
168    // Initialization
169    //
170    
171    if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
172    
173    switch (fCandidateType) {
174    case 0: 
175       {
176          AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
177          copyfCutsDzero->SetName("AnalysisCutsDzero");
178          PostData(2, copyfCutsDzero);  // Post the data
179       } break;
180    case 1: 
181       {
182          AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
183          copyfCutsDstar->SetName("AnalysisCutsDStar");
184          PostData(2, copyfCutsDstar); // Post the cuts
185       } break;
186    default: return;
187    }
188    
189    return;
190 }
191
192 //_______________________________________________________________________________
193
194 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
195
196    //
197    // output 
198    //
199    
200    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
201    
202    fOutput = new TList(); fOutput->SetOwner();
203    DefineHistoForAnalysis(); // define histograms
204    
205    if (fCandidateType==kD0toKpi){
206       fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
207       fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0); 
208    }
209    
210    if (fCandidateType==kDstartoKpipi) {
211       fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
212       fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); 
213    }
214    
215    fCandidateArray->SetOwner();
216    fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
217    
218    //this is used for the DStar side bands and MC!
219    fSideBandArray->SetOwner();
220    fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
221   
222    PostData(1, fOutput);
223    PostData(3, fCandidateArray);
224    PostData(4, fSideBandArray);
225  
226    return;
227 }
228
229 //_______________________________________________________________________________
230
231 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
232    //
233    // user exec
234    //
235    // Load the event
236    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
237    
238    TClonesArray *arrayDStartoD0pi = 0;
239    if (!aodEvent && AODEvent() && IsStandardAOD()) {
240       
241       // In case there is an AOD handler writing a standard AOD, use the AOD 
242       // event in memory rather than the input (ESD) event.    
243       aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
244       
245       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
246       // have to taken from the AOD event hold by the AliAODExtension
247       AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
248       if(aodHandler->GetExtensions()) {
249          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
250          AliAODEvent *aodFromExt = ext->GetAOD();
251          arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
252       }
253    } else {
254       arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
255    }
256    
257    if (!arrayDStartoD0pi) {
258       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
259       return;
260    } else {
261       AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
262    }
263    
264    TClonesArray* mcArray = 0x0;
265    if (fUseMCInfo) {
266       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
267       if (!mcArray) {
268          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
269          return;
270       }
271    }
272    
273    //Histograms
274    TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
275    TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
276    TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
277    TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
278    
279    TH1F* hPtPion=0x0;
280    if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
281    hstat->Fill(0);
282    
283    // fix for temporary bug in ESDfilter 
284    // the AODs with null vertex pointer didn't pass the PhysSel
285    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
286    
287    //Event selection
288    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
289    //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
290    if (!iseventselected) return;
291    hstat->Fill(1);
292    
293    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
294    if(fUseReco) hstat->Fill(2, nD);
295    
296    //D* and D0 prongs needed to MatchToMC method
297    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
298    Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
299    
300    //D0 from D0 bar
301    Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K 
302    Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi 
303    
304    Int_t iCand =0;
305    Int_t iSBCand=0;
306    Int_t isSelected = 0;
307    AliAODRecoDecayHF *charmCand = 0;
308    AliAODRecoCascadeHF *dstar = 0;
309    AliAODMCParticle *charmPart = 0;
310    Bool_t isMCBkg=kFALSE;
311    
312    Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
313    
314    Int_t mcLabel = -9999;
315    Int_t pdgMeson = 413;
316    if (fCandidateType==kD0toKpi) pdgMeson = 421;
317    
318    //clear the TClonesArray from the previous event
319    fCandidateArray->Clear();
320    fSideBandArray->Clear();
321    
322    for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
323       charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
324       if (!charmCand) continue;
325       
326       if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
327       
328       if (fUseMCInfo) { // Look in MC, try to simulate the z
329          if (fCandidateType==kDstartoKpipi) {
330             mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
331          }
332          
333          if (fCandidateType==kD0toKpi) 
334             mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
335          
336          if (mcLabel<=0) isMCBkg=kTRUE;
337          else hstat->Fill(2);
338          if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
339       }
340       
341       Double_t ptD = charmCand->Pt();
342       
343       // region of interest + cuts
344       if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;    
345       
346       if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
347          //select by track cuts the side band candidates (don't want mass cut)
348          isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent); 
349          if (!isSelected) continue;
350          //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
351          Int_t bin = fCuts->PtBin(ptD);
352          if(bin<0 || bin>=fCuts->GetNPtBins()) {
353             AliError(Form("Pt %.3f out of bounds",ptD));
354             continue;
355          }
356          //if data and Dstar from D0 side band
357          if (((dstar->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((dstar->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){       
358             
359             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
360             iSBCand++;
361          }
362       }
363       //candidate selected by cuts and PID
364       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
365       if (!isSelected) continue;
366       
367       //for data and MC signal fill fCandidateArray
368       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
369          // for data or MC with the requirement fUseReco fill with candidates
370          if(fUseReco) {
371             if (fCandidateType==kDstartoKpipi){
372                new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
373                AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
374             } else{
375                new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
376                //Printf("Filling reco");
377             }               
378             
379             hstat->Fill(3);
380          }
381          // for MC with requirement particle level fill with AliAODMCParticle
382          else if (fUseMCInfo) {
383             new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
384             //Printf("Filling gen");
385             hstat->Fill(3);
386          }
387          
388          iCand++;       
389       }
390       //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
391       else if(fUseReco){
392          if (fCandidateType==kDstartoKpipi){
393             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
394          }
395          if (fCandidateType==kD0toKpi){
396             new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
397          }
398          iSBCand++;
399       }
400       
401       
402       Double_t masses[2];
403       if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
404          
405          //softpion from D* decay
406          AliAODTrack *track2 = (AliAODTrack*)dstar->GetBachelor();  
407          
408          // select D* in the D0 window.
409          // In the cut object window is loose to allow for side bands
410          
411          
412          // retrieve the corresponding pt bin for the candidate
413          // and set the expected D0 width (x3)
414          //    static const Int_t n = fCuts->GetNPtBins();
415          Int_t bin = fCuts->PtBin(ptD);
416          if(bin<0 || bin>=fCuts->GetNPtBins()) {
417             AliError(Form("Pt %.3f out of bounds",ptD));
418             continue;
419          }
420          
421          AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
422          //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
423          if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {      
424             masses[0] = dstar->DeltaInvMass(); //D*
425             masses[1] = 0.; //dummy for D*
426             
427             //D*  delta mass
428             hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
429             
430             // D* pt and soft pion pt for good candidates               
431             Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
432             Double_t invmassDelta = dstar->DeltaInvMass();
433             if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
434          }
435       } //Dstar specific
436       
437       if (fCandidateType==kD0toKpi) { //D0->Kpi
438          
439          //needed quantities
440          masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
441          masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
442          hstat->Fill(3);
443          
444          // mass vs pt
445          if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
446          if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
447       } //D0 specific
448       
449       charmCand = 0;
450       isMCBkg=kFALSE;
451    } // end of D cand loop
452    
453    hnCandEv->Fill(fCandidateArray->GetEntriesFast());
454    if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
455       Int_t nsbcand=fSideBandArray->GetEntriesFast();
456       hstat->Fill(4,nsbcand);
457       hnSBCandEv->Fill(nsbcand);
458    }
459    Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
460    PostData(1, fOutput);
461    PostData(3, fCandidateArray);
462    PostData(4, fSideBandArray);
463   
464    return;
465 }
466
467 //_______________________________________________________________________________
468
469 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
470 {
471    // The Terminate() function is the last function to be called during
472    // a query. It always runs on the client, it can be used to present
473    // the results graphically or save the results to file.
474    
475    Info("Terminate"," terminate");
476    AliAnalysisTaskSE::Terminate();
477    
478    fOutput = dynamic_cast<TList*>(GetOutputData(1));
479    if (!fOutput) {
480       printf("ERROR: fOutput not available\n");
481       return;
482    }
483    
484    return;
485 }
486
487 //_______________________________________________________________________________
488
489 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
490 {
491    //
492    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
493    //
494    
495    Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
496    
497    // compute the Delta mass for the D*
498    if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
499    
500    
501    fMinMass = mass - range;
502    fMaxMass = mass + range;
503    
504    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
505    if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
506    
507    return;
508 }
509
510 //_______________________________________________________________________________
511
512 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
513 {
514    //
515    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
516    //
517    
518    if (uplimit>lowlimit) {
519       fMinMass = lowlimit;
520       fMaxMass = uplimit;
521    } else {
522       printf("Error! Lower limit larger than upper limit!\n");
523       fMinMass = uplimit - uplimit*0.2;
524       fMaxMass = uplimit;
525    }
526    
527    return;
528 }
529
530 //_______________________________________________________________________________
531
532 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
533 {
534    //
535    // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
536    //
537    
538    if (nptbins>30) {
539       AliInfo("Maximum number of bins allowed is 30!");
540       return kFALSE;
541    }
542    
543    if (!width) return kFALSE;
544    for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
545    
546    return kTRUE;
547 }
548
549 //_______________________________________________________________________________
550
551 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
552 {
553    //
554    // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
555    //
556    
557    // Statistics 
558    TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
559    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
560    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
561    if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
562    else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
563    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
564    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");  
565    hstat->SetNdivisions(1);
566    fOutput->Add(hstat);
567    
568    TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
569    fOutput->Add(hnCandEv);
570    
571    // Invariant mass related histograms
572    const Int_t nbinsmass = 200;
573    TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
574    hInvMass->SetStats(kTRUE);
575    hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
576    hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
577    fOutput->Add(hInvMass);
578    
579    if (fCandidateType==kDstartoKpipi) {
580       TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
581       fOutput->Add(hnSBCandEv);
582       
583       TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
584       hPtPion->SetStats(kTRUE);
585       hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
586       hPtPion->GetYaxis()->SetTitle("entries");
587       fOutput->Add(hPtPion);
588    }
589    
590    return kTRUE; 
591 }